Method and controller for controlling an aircraft by improved direct lift control

ABSTRACT

The present invention discloses a method and flight controller for controlling an aircraft, comprising the following step: —determining the pitch control input δE for at least one pitch moment generator element, based on the lift control input δF for at least one lift generator element; wherein the step of determining the pitch control input δE based on the lift control input δF includes the step of determining a feed forward filter output Fq,δE,δF(q).

CROSS-REFERENCE TO RELATED APPLICATION DATA

This application is a US National Stage Application ofPCT/EP2018/064932, filed Jun. 6, 2018, titled Method and Controller forControlling an Aircraft by Improved Direct Lift Control, which claimsthe benefit of and priority to Austrian Application No. A 60048/2017,filed Jun. 7, 2017, the disclosures of which are incorporated herein intheir entirety.

FIELD OF THE INVENTION

The present invention relates to a method and controller for controllingan aircraft, particularly for controlling the vertical acceleration. Theinvention also relates to an aircraft and a simulator using theinventive method or controller.

Particularly the invention relates to a so called direct lift control.To this end, lift generating elements like flaps, slats, spoilers,micro-thruster or the like are actuated to generate the so called directlift. The principles of the direct lift concept are known to the manskilled in the art.

A controller of an aircraft may use the so-called path following control(PFC) concept for a fixed wing aircraft. Prior concepts of direct liftcontrol suffer from undesired resonances due to angle of attackoscillations. The present invention proposes a control concept for thevertical acceleration allowing to more than double a controller'sbandwidth compared to classical approaches. By using a dynamic flapactuation, the lift, and therefore vertical acceleration, can bemanipulated more quickly than by rotating the whole aircraft by means ofthe elevator. Thus, in the disclosed concept, flap deflection is used asprimary control input, and the elevator serves for maintaining the wingsin their operational range and returning the flaps to their neutralposition.

The disadvantages of classical lift generation by means of elevatordeflection are caused by the necessary rotation of the aircraft. Therotation of the aircraft to change the angle of attack and therebyindirectly increase the lift causes a delay of lift generation whichlimits the possible bandwidth of lift control. Furthermore, due to therotation the pitch angle of the aircraft varies for indirect liftgeneration. In consequence, angular accelerations lead to undesiredaccelerations at cabin areas which are located distant from the centerof gravity. In this context, direct lift control aims to overcome thesedrawbacks by directly generating lift without requiring the rotation ofthe aircraft. Unfortunately, for state-of-the-art direct lift control anangle of attack oscillation appears in the dynamic region of theresonance frequency of the short-period mode. This angle of attackvariation is not caused by additional pitching moments, which may begenerated by direct lift control surfaces, but due to the aircraft'smotion caused by the lift variation. The resulting angle of attackvariation generates a lift variation which opposes the intended directlift. Therefore, an antiresonance occurs which lowers the effectivenessof state-of-the-art direct lift control.

PRIOR ART

U.S. Pat. No. 4,261,537 and US 2016/023749 A1 disclose a prior artdirect lift control.

US 2009/157239 A1 discloses an elevator gain shaping providing an inputsignal for the elevator control.

SUMMARY OF THE INVENTION

It is an object of the invention to provide a method and controller forcontrolling an aircraft or an aircraft simulator by use of a direct liftcontrol.

The object of the present invention is solved by a method according toclaim 1 a computer program product according to claim 9, and a flightcontroller according to claim 10. The depending claims are directed toembodiments of the invention.

The disclosed method for controlling an aircraft includes the step ofdetermining the pitch control input δ_(E) for at least one pitch momentgenerator element, based on the lift control input δ_(F) for at leastone lift generator element, wherein the step of determining the pitchcontrol input δ_(E) based on the lift control input δ_(F) includes thestep of determining a feed forward filter output F_(q,δE,δF)(q) of thelift control input or state of lift generating element.

The disclosed method and flight controller concern an aircraft, whichadvantageously includes at least one lift generator element, the directeffect of which on the lift of the aircraft is faster than that of theelevator or other types of pitch moment generator elements. If the liftgenerator element can alter the lift only to one direction, i.e. onlyincrease or only decrease lift, a partial pre-actuation of the liftgenerator element allows subsequent alteration of lift in bothdirections. E.g. if spoiler deflection only allows for lift decrease, apre-deflection of the spoiler allows subsequent deflection of thespoiler in both directions and thereby alter the lift in bothdirections. The lift control input δ_(F) can be commanded by the pilot,a flight controller, an autopilot or the like.

In a first embodiment, at least one lift generator element is anaerodynamic surface of the aircraft adapted to modify the lift of theaircraft. The aerodynamic surface of the aircraft is preferably one ofthe following elements: a spoiler of the aircraft; a flaperon of theaircraft; a flap of the aircraft; a slat of the aircraft.

In a second embodiment, at least one lift generator element is amicro-thruster or propulsion unit adapted to generate an actionmodifying the lift of the aircraft or directly generating a force inlift direction.

In a third embodiment, at least one lift generator element is a wing,which is capable of changing the incidence angle and/or the wing shape.

For the aforementioned embodiments the pitch moment generator element ispreferably one of the following: an elevator of the aircraft; ahorizontal stabilizer of the aircraft; a canard foreplane of theaircraft; a micro-thruster of the aircraft; a propulsion unit of theaircraft; a wing, which is capable of changing the incidence angleand/or the wing shape.

The invention also discloses a method for controlling an aircraft,comprising the steps of determining the commanded elevator angle δ_(E)based on the flap element angle δ_(F), wherein the step of determiningthe commanded elevator angle δ_(E) based on the flap element angle δ,includes the step of determining a feed forward elevator deflectionF_(q,δE,δF)(q) of the elevator angle. The method may control thevertical acceleration of the aircraft. The flap element angle may be thecommanded flap element angle and/or the actual flap element angle. Inthis embodiment, the flap element angle δ, is the lift control input andthe commanded elevator angle δ_(E) is the pitch control input. The flapelement angle δ_(F) can be commanded by the pilot, a flight controller,an autopilot or the like.

The invention also discloses a flight controller adapted to control atleast one of an aircraft and a component of a flight simulator. Theflight controller comprises a feed forward filter adapted to determinethe elevator angle δ_(E) based on the flap element angle δ_(F), whereinthe feed forward filter F_(q,δE,δF)(q) determines a feed forwardelevator deflection based on the flap element angle δ_(F).

The benefit of including the feed forward filter is that the angle ofattack oscillation, which is caused by direct lift generation and leadsto the antiresonance, is reduced and consequently the effectiveness ofthe direct lift control is improved. Additionally, using the flapelements as primary control input further raises the achievablebandwidth of the disclosed control method. The improved dynamic responseof the system amongst others facilitates higher flight path precision,disturbance rejection, turbulence alleviation and touchdown precisionduring landings.

The embodiments are described below with reference to the method andflight controller for the sake of brevity.

The flap element angle may be a flap angle, a flaperon angle, a spoilerangle, a slat angle or the like. The control concept will be explainedand demonstrated for an aircraft comprising two flaps and two flaperons.

Determining the elevator angle δ_(E) based on the flap element angle δ,may comprise a high pass filtering the flap element angle δ, by a highpass, a high pass filtering the flap element angle δ_(F) by alead-filter or a high pass filtering the flap element angle δ_(F) by ahigh pass of second order.

In one embodiment the feed forward elevator deflection F_(q,δE,δF)(q) isdetermined by the following formula:

${F_{q,{\delta\; E},{\delta\; F}}(q)} = \frac{{b_{1}q} + {b_{2}q^{2}}}{q^{2} + {2\;\xi_{F}\omega_{F}q} + \omega_{F}^{2}}$whereinq is the Laplace variable;b₁ is a first optimization parameter;b₂ is a second optimization parameter;ξ_(F) is a third optimization parameter; andω_(F) is a fourth optimization parameter.ξ_(F) may be a filter damping factor. ω_(F) may be a resonancefrequency.

The parameter b₁, b₂, ξ_(F) and ω_(F) may be optimized by minimizing thefollowing cost function:J=Σ _(ω) _(k) _(∈Ω) |{tilde over (P)} _(q,δ) _(F) (jω _(K))−{tilde over(P)} _(q,δ) _(F) ^(des)(jω _(k))|wherein {tilde over (P)}_(δ) _(F) is an optimal fit to a desired plant{tilde over (P)}_(δ) _(F) ^(des) by tuning the parameter vector [b₁ b₂ξ_(F) ω_(F)].

Determining the commanded elevator angle δ_(E) based on the candidateelevator angle δ_(E) ^(C) and on the flap element angle δ_(F) mayinclude adding the feed forward elevator deflection F_(q,δE,δF)(q) to acandidate elevator angle δ_(E) ^(C). The flight controller may comprisean adder for adding the feed forward elevator deflection F_(q,δE,δF)(q)to the candidate elevator angle δ_(E) ^(C). The candidate elevator angleδ_(E) ^(C) can be commanded by the pilot, a flight controller, anautopilot or the like.

The method and flight controller may determine the commanded flapelement angle δ_(F) ^(des) and determine the candidate elevator angleδ_(E) ^(C) based on the flap element angle δ_(F) and the desired flapelement angle δ_(F) ^(des) of the aircraft, such as by a firstcontroller. The desired flap element angle δ_(F) ^(des) may be commandedby a pilot or any flight computer.

Determining a candidate elevator angle δ_(E) ^(C) may include a PIcontrol based on the flap element angle δ_(F) and the desired flapelement angle δ_(F) ^(des) of the aircraft. The first controller may bea PI controller.

The method and flight controller may comprise or execute the followingsteps, e.g. by a second controller: determining the actual verticalacceleration a_(zB) of the aircraft, determining the desired verticalacceleration a_(zB) ^(des) of the aircraft, and controlling the flapelement angle δ_(F) based on the actual vertical acceleration a_(zB) ofthe aircraft and desired vertical acceleration a_(zB) ^(des) of theaircraft. The method may implement a PI control. The second controllermay be a PI controller. The desired desired vertical acceleration a_(zB)^(des) may be commanded by a pilot or any flight computer.

The control based on the commanded flap element angle δ_(F) and thedesired flap element angle δ_(F) ^(des) of the aircraft has a lowerbandwidth than the control based on the actual vertical accelerationa_(zB) of the aircraft and desired vertical acceleration a_(zB) ^(des)of the aircraft. The first controller may have a lower bandwidth thanthe second controller.

The vertical acceleration a_(zB) of the aircraft and the desiredvertical acceleration a_(zB) of the aircraft are normalized by thefollowing formula, such as by a normalization unit:

$c_{zB} = {\frac{a_{aB}}{V_{A}^{2}}\frac{V_{ref}^{2}}{a_{ref}}}$whereinc_(zB) is the coefficient of normalized acceleration;V_(A) is the actual air speed;V_(ref) is the reference air speed; anda_(ref) is the vertical reference acceleration.

In one embodiment, the actual air speed V_(A) and the actual verticalacceleration a_(zB) may be provided by a sensor of an aircraft.

In another embodiment, the coefficient of normalized accelerationc_(zB), the actual air speed V_(A), and the actual vertical accelerationa_(zB) may be provided by a simulator.

The invention also discloses a computer program product that when loadedinto a memory of a computer having a processor executes the abovemethod. The invention also discloses an aircraft comprising the aboveflight controller. The invention also discloses a simulator comprisingthe above flight controller.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention is now described in detail with reference to theaccompanying drawings showing exemplary and not restricting embodimentsof the invention, wherein:

FIG. 1 shows a schematic top view of an aircraft implementing theinvention;

FIG. 2 shows a table for an equilibrium state for P_(ω) _(yB)-identification;

FIG. 3 shows a table for an equilibrium state for P_(ω) _(xB)-identification;

FIG. 4 shows a table for an equilibrium state for P_(a) _(yB)-identification;

FIG. 5 shows Bode plots of P_(q,δ) _(E) and P_(q,δ) _(F) from δ_(E) andδ_(F), respectively, to the response −c_(zB).

FIG. 6 shows a step response of −c_(zB) due to step input of δ_(E) andδ_(F), respectively;

FIG. 7 depicts a trajectory for sinusoidal acceleration a_(zB) andcancellation effect due to short-period mode lag;

FIG. 8 illustrates an original plant P_(q,δ) _(F) , desired plant {tildeover (P)}_(q,δ) _(F) ^(des), filter F_(q,δ) _(E) _(,δ) _(F) , andresulting plant {tilde over (P)}_(q,δ) _(F) ;

FIG. 9 shows an equilibrium state for P_(δ) _(F) and {tilde over(P)}_(δ) _(F) -identification;

FIG. 10 depicts measured response −c_(zB) to the chirp of δ_(F) withoutand with F_(δ) _(E) _(,δ) _(F) resulting in {tilde over (P)}_(δ) _(F) ;

FIG. 11 shows a closed loop step response of −c_(zB) and δ_(F) for−c_(zB) ^(des)=1;

FIG. 12 illustrates a structure of a c_(zB)-control including C_(c)_(zB) , C_(δ) _(F) and a designed filter F_(δ) _(E) _(,δ) _(F) ;

FIG. 13 depicts a table showing an equilibrium state for P_(δ) _(F)_(,δ) _(E) -identification;

FIG. 14 illustrates a response of δ_(E) and δ_(F) to an outputdisturbance δ_(F)=1;

FIG. 15 depicts a setpoint change from c_(zB,0) to c_(zB,1); and

FIG. 16 shows a table for an equilibrium state for P_(a) _(x)_(B)-identification.

DETAILED DESCRIPTION OF THE DRAWINGS

Reference is made to FIG. 1. A fixed wing aircraft 100 used forvalidation of the concept is a model of a Piper PA18 Super Cub with aspan of 1300 mm. The original version has a span of 10.7 m, thus, thescaling factor is 1:8.25. This type of airplane is known for very goodtake-off and landing performance. Therefore, it is also considered to bea “Short Take-Off and Landing” (STOL) aircraft.

The fixed wing aircraft 100 comprises wings 102 and a tail 104 having anelevator 106 and a rudder 108. The wings 102 comprise ailerons 110, 111and flaps 114. A motor 116 drives a propeller 118. The aircraft 100 iscontrolled by a flight controller 120, which is connected to the firstservo 124 actuating the left aileron 110 and to the second servo 130actuating the elevator 106. The flight controller 120 is furtherconnected to the third servo 138 controlling the motor 116. A forthservo 128 actuating the rudder 108 is connected to the flight computer120. The fifth servo 126 actuating the flaps 114 is also connected tothe flight controller 120. The sixth servo 136 actuating the rightaileron 111 is connected to the flight controller 120. The aircraftcomprises a Pitot-static-tube 122 for measuring the air speed, whereinthe Pitot-static-tube 122 is operatively coupled with the flightcontroller 120. The aircraft 100 further comprises a GNSS and magneticflux sensor 134 operatively coupled with the flight controller 120 fordetermining the absolute position of the aircraft 100.

A special feature of this particular model aircraft are vortexgenerators 103, which are located in the first quarter of the uppersurface of the wings 102. By purposeful disturbance of the boundarylayer, an early transition from laminar airflow to turbulent airflow istriggered. This entails an intensified exchange of lower energetic airof the boundary layer with higher energetic air of the outside air flow.

Therefore, the stall characteristic of the wing can be improved,particularly the boundary layer stays attached for higher angles ofattack. As rapid flap deflections will be used, the vortex generatorsmay improve the response to intended lift manipulations. The interfacebetween the physical world and the mathematical quantities, which areinputs and outputs of the calculations for the flight control, is givenby sensors and actuators.

As to the air speed sensor 122, by means of a pitot-static tube, thedifferential pressure between total pressure p_(t) and static pressurep_(s), called dynamic pressure

${q = {{p_{t} - p_{s}} = {\frac{\rho}{2}V_{A}^{2}}}},$is determined, with the air density ρ, and the airspeed V_(A). Thus, theairspeed can be calculated as

$V_{A} = {\sqrt{\frac{2\left( {p_{t} - p_{s}} \right)}{\rho}}.}$The root-characteristic of the measurement principle implies that forsmall V_(A) and therefore small pressure difference q=p_(t)−p_(s), smallvariations of the measured pressure Δq yield large airspeed variationsΔV_(A). Therefore, noise of the pressure measurement has got high impacton the measured airspeed on the ground, i. e., for V_(A)≈0 m/s. However,in-flight, for higher airspeeds, the determination of V_(A) by means ofa pressure measurement becomes more reliable.

The purpose of measuring the airspeed V_(A) is to have information aboutthe aerodynamic effects, which depend on the speed in the air ratherthan the speed over ground. For an accurate measurement, thepitot-static tube must be placed where the air flow is ideallyundisturbed. The pitot-static tube has to be placed outside the boundarylayer, i. e., not too close to the aircraft's surface, where the airflow slows down. The speed of the airflow varies on the wing's upper andlower surface as well as at different fuselage positions, depending onthe curvature of the surfaces and how the airflow approaches. Thus, aposition without significant air stream deflection is endeavored.Considering a single motor aircraft like the Piper PA18, the airflow atthe fuselage is highly affected by the propeller's slipstream.Therefore, an airspeed measurement close to the fuselage would bedisturbed by the thrust setting. Placing the pitot-static tube too farto the wing tip, e. g., at a distance d_(tip), would cause coupling ofangular rates about the vertical body axis ω_(zB), due to the superposedtangential velocity ΔV=d_(tip)ω_(zB).

Bearing in mind these criteria, the pitot-static tube is installedprotruding from the boundary layer, in front of the wing beforesignificant airflow deflection, outside the propeller's slipstream, butstill not too far towards the wing tip.

As to the flight controller 120 comprising the accelerometer, gyroscopesand altimeter 136 by means of MEMS technology, multiple accelerometers136 measuring the acceleration along the body axes are installed. Thedifference of acceleration between the body to be measured, i. e., theaircraft 100, and a reference mass, which is part of the MEMS-structure,is determined. Thus, MEMS accelerometers are not capable of measuringgravity, which is a body force, and therefore, acts the same way on bothbodies. The measurement principle of the gyroscopes 136 is based on theintentionally excited oscillation of a MEMS-structure.

An angular rate about a specific axis causes a secondary oscillation dueto Coriolis effects, which is measured. By measuring in three orthogonaldirections, the accelerations and angular rates of all three body axesare measured. Finally, the barometric altitude is measured by abarometer 136, which measures the static pressure. Due to anapproximately exponential pressure drop as a function of altitude, thebarometric altitude can be calculated from the pressure measurement.

The Pixhawk Flight-Controller 120, which comprises the above sensors136, is placed inside the fuselage 132, close to the center of gravity.Besides the available space at this position, placing the accelerometers136 at the center of gravity provides the advantage of not measuringparasitic accelerations due to angular movements.

The combined module 134 delivers GPS (GNSS) position data as well as 3Dmagnetic flux information. The GPS position is determined by theposition and time information received from at least four satellites. Bydistance calculation based on the time of flight of the received signalsand trilateration, the desired information of position and time isobtained. By means of anisotropic magneto-resistive sensors 134, the3-axes digital compass converts any incident magnetic field in thesensitive axes directions to a differential voltage output whichrepresents the measured magnetic flux. For best reception of the GPSsignal, placing the module on the top face of the aircraft isendeavored. The magnetic flux should be placed at a spot of low magneticinterference, caused by the BLDC-motor 116 and motor controller 138,which are positioned in the front part of the aircraft. Therefore, thesensor unit 134 is placed at the rear part of the top face.

The actuators of the model aircraft are mainly the servos 124, 126, 128,130 136, 138, which deflect the control surfaces, and the throttle inputfor the BLDC motor controller. All control surfaces are configured tohave a maximum deflection of about 25°. The input signals for elevatordeflection δ_(E), and rudder deflection δ_(R), are directly used as theservo inputs for the elevator servo 130 and rudder servo 128,respectively. For positive δ_(E), the elevator moves up, for positiveδ_(R) the rudder moves to the right. The input signal for ailerondeflection δ_(A) actuates the servos 124, 136 of the left and rightaileron in opposite directions. For positive δ_(A), the right aileron111 moves up and the left aileron 110 moves down. The input signal forflap deflection δ_(F) also actuates the aileron servos, and additionallythe servo of the flaps. For positive δ_(F) all three servos move down tothe same extent. Thus, the ailerons are used as so called “flaperons”,i. e., control surfaces, which are used as both, aileron and flaps.Therefore, the input for flap deflection 6F influences the curvature ofthe whole wing, instead of only the middle part. This results in twofavorable effects. Firstly, the achieved increase or decrease of lift isgreater. Secondly, the parasitic effect of a pitching moment is reduced,which shall be briefly outlined. The control surface deflectionincreases the lift in the rear part of the wing, which is behind thecenter of gravity. Thus, in principle a pitch down moment is generatedby flap deflection. However, for the inner sector, where the horizontalstabilizer is located aft of the wing, another effect dominates, whichcauses a pitch up moment. Due to the increased lift, the flap deflectiongenerates increased downwash, i. e., a downwards movement of theairflow. Subsequently, a downwards force at the elevator is generated,which generates a pitch up moment. As a consequence, flap deflection atthe outer wing causes a pitch down moment, while at the inner wing thepitch up moment dominates. By using both, the outer and the innercontrol surfaces, large parts of the effects cancel out.

The flight controller Pixhawk 120 is an open hardware and open softwareproject, which is the outcome of a development for many years. Thestandard control algorithms of the flight stack ArduPilot aren't used atall. However, the rest of the framework, like scheduling, sensor andactuator communication, and the flight state estimations are usedwithout further modifications.

The development of the control concept can be performed with theassistance of various simulation tools. The advantage of simulation isto have full access to all quantities, the easy adjustment ofparameters, and the possibility to step by step add effects such asnoise and disturbances, e. g., wind. In addition to engineering toolslike Matlab/Simulink, the following tools can be used. By use ofSoftware-In-The-Loop (SITL) simulation with X-Plane not only theaircraft dynamics can be simulated, but also the behavior of the flightcontroller Pixhawk. Thus, instead of Matlab functions, a C++ code can beimplemented, which is already part of the flight stack ArduPilot. Byconnecting the physical remote control to the computer, which runs theSITL and X-Plane, the actual flight conditions of a real test flight canbe simulated.

The last step to validate the function of the control concept beforereal test flights is to include also the real hardware(Hardware-In-The-Loop (HITL) simulation with X-Plane). Therefore, viaUSB the flight controller Pixhawk can be connected to the computer whichruns X-Plane. The control inputs and flight states are exchanged via UDPpackages. Instead of values of the real sensors, e. g., the calculatedposition of the GPS module, the values from simulation are used.

The simulation delivers good results which are validated with the realset-up. The basic control concept simulated with a real size Cessna C172is the same as for the test flights with the model aircraft Piper PA18.The inner loop controllers and the geometric properties are adjusted tofit to the real size airplane. As the PFC is based on abstract kinematicquantities such as accelerations, the concept can be appliedindependently of the aircraft size. Only the achievable bandwidth of theinner loop control and characteristics like stall speed, maximumairspeed, and acceleration restrictions depend on the specific type ofaircraft. Simulations with an Airbus A320 and a Boeing B737 show thatthe control concept can also be applied on airliner effectively.

For modelling the fixed wing aircraft dynamics, the aircraft isconsidered to be a rigid body. Two reference frames are used. Firstly,the inertial frame (Index I) is defined according to the North-East-Downcoordinate system, which is also known as local tangent plane (LTP). Ithas a fixed origin on the surface of the earth with the x-axis to theNorth, y-axis to the East, and the z-axis according to a right-handedsystem pointing Down. Secondly, the body reference frame (Index B) ischosen to have its origin at the center of gravity of the aircraft withthe x-axis in longitudinal direction, the y-axis in lateral directionpointing to the right wing, and the z-axis in vertical directionpointing down.

By representing the acting forces f_(B) and torques τ_(B) in the bodyreference frame, the state space model of a rigid body can be found as

$\begin{matrix}{{\overset{.}{r}}_{I} = {{R_{I}^{B}\left( {\varphi,\theta,\psi} \right)}v_{B}}} & (1) \\{{\overset{.}{v}}_{B} = {\frac{1}{m}\left( {f_{B} - {\omega_{B} \times \left( {mv}_{B} \right)}} \right)}} & (2) \\{{\overset{.}{\vartheta}}_{RPY} = {{R_{RPY}^{B}\left( {\varphi,\theta,\psi} \right)}\omega_{B}}} & (3) \\{{{\overset{.}{\omega}}_{B} = {I_{B}^{- 1}\left( {\tau_{B} - {\omega_{B} \times \left( {I_{B}\omega_{B}} \right)}} \right)}},} & (4)\end{matrix}$where r_(I)=[x_(I) y_(I) z_(I)]^(T) is the position represented in theinertial reference frame, v_(B) the velocity represented in the bodyreference frame, ϑ_(RPY)=[φ θ ψ]^(T) the orientation expressed inRoll-Pitch-Yaw representation ω_(B)=[ω_(xB) ω_(yB) ω_(zB)]^(T)=[p qr]^(T) the angular rates represented in the body reference frame, R_(I)^(B) the transformation matrix from the body frame to the inertialframe, R_(RPY) ^(B) the transformation matrix from the body angularrates to Roll-Pitch-Yaw rates, I_(B) the moments of inertia representedin the body reference frame, and m the mass of the aircraft.

The body forces f_(B) basically consist of lift L, drag D, side forceSF, and gravity, which need to be oriented correctly. One example, whichwill be used below is the equation for modelling the lift force L

$\begin{matrix}{{L = {c_{L}\frac{\rho}{2}V_{A}^{2}S}};} & (5) \\{{c_{L} \approx {c_{L\; 0} + {c_{L,\alpha}\alpha} + {c_{L,q}q} + {c_{L,\delta_{E}}\delta_{E}} + {c_{L,\delta_{F}}\delta_{F}}}},} & (6)\end{matrix}$where ρ denotes the air density, V_(A) the true airspeed, S the wingarea, and c_(L) the lift coefficient. The latter is usually approximatedusing an affine representation of the angle of attack α, the angularrate q, and the deflections of the elevator δ_(E) and of the flapsδ_(F).

In general, the dynamics of a specific aircraft design result from thedependencies of the acting forces f_(B) and torques τ_(B) on the 12rigid body states r_(I), v_(B), ϑ_(RPY), ω_(B), and

the physical inputs δ_(A), δ_(E), δ_(F), δ_(R), δ_(T) being given byaileron deflection, elevator deflection, flap deflection, rudderdeflection and throttle actuation. Aileron deflection δ_(A) in firstplace results in a roll moment, a torque about the body x-axis. Thus, itis used to influence the roll rate and subsequently the roll angle.Parasitic effects are drag variation and yawing moments.

Elevator deflection δ_(E) in first place results in a pitch moment, atorque about the body y-axis. Thus, it is used to influence the angle ofattack and subsequently the lift, or the pitch rate and subsequently thepitch angle. Parasitic effects are forces in direction of the verticalaxis, which are necessary to generate the pitch moment and add to thelift of the wings. Depending on whether the horizontal stabilizer is atailplane or a foreplane, the parasitic force decreases or increases thelift. In the case of a frequently used tailplane, the resulting liftforce results to have a non minimum-phase characteristic.

Flap deflection δ_(F) in first place results in additional lift and tosome extent additional drag. Parasitic effects are pitch moments. Thusconventionally, flaps are used during take-off and landing to increaselift for reducing the airspeed. Furthermore, the additional drag can bedesirable to decelerate and it allows for steeper approach angleswithout exceeding speed limitations. In this disclosure, flap deflectionis used as primary control input for the variation of lift and thereforevertical acceleration.

Rudder deflection δ_(R) in first place results in a yaw moment, a torqueabout the body z-axis. Thus, it is used to influence the side slip angleand subsequently the side force, or the yaw rate and subsequently theyaw angle. Parasitic effects are roll moments and a parasitic sideforce, which is necessary to obtain the yaw moment.

In cruise flight the rudder is primarily used to maintain a coordinatedflight, i. e., zero side force, indicated in the cockpit by aninclinometer called slip indicator. In the final phase of a landingprocedure, the rudder is used for yaw control to align the aircraft withthe runway centerline.

Throttle actuation δ_(T) in first place results in a thrust variation,i. e., force in longitudinal direction. Parasitic effects are roll,pitch, and yaw moments and lift variation.

The inputs δ_(A), δ_(E), δ_(F), δ_(R) are normalized to [−1, 1]. δ_(T)is normalized to [0, 1].

The aircraft dynamics of order twelve can be separated into 6longitudinal states, according to dynamic modes in the longitudinalplane, and six lateral-directional states, resulting in dynamic modesorthogonal to the longitudinal plane.

The phugoid mode (longitudinal, order 2) is basically an interchange ofkinetic energy and potential energy causing airspeed and altitudevariations. It can be illustrated by observing one oscillation period.Starting when the altitude is highest, the velocity is smallest, whichcauses a lack of lift. Thus, altitude decreases, which causes a raise invelocity and therefore lift. In further consequence, due to theincreasing lift the aircraft starts to climb again. Finally, thealtitude is maximal at the end of one period again. The pitch angle alsochanges due to the up- and downwards varying flight path, while theangle of attack variation is only small. The period length is in theorder of 60 s for civil airplanes and 10 s for smaller UAVs, mainlydepending on the airspeed, and usually the phugoid mode is weaklydamped.

The short-period mode (longitudinal, order 2) is essentially an angle ofattack oscillation caused by a resetting torque towards the equilibriumangle of attack, which can be influenced by elevator deflection δ_(E).The period length is in the order of a few seconds for civil airplanesand of 0.5 s for smaller UAVs. This mode is usually considerably damped,i. e., it shows small overshoots only.

The roll mode (lateral, order 1) is caused by the damping of the rollrate, the angular rate about the longitudinal axis. For constant ailerondeflection δ_(A), the roll rate is approaching a magnitude where inputtorque caused by aileron deflection equals the counteracting dampingtorque.

The spiral mode (lateral, order 1) is describing the behavior of theaircraft to regain straight flight conditions after a roll angleperturbation or to further increase the roll angle, resulting in aspiral dive. Depending on the aircraft design this mode can be stable orunstable. However, it exhibits a long time constant for usual aircraftdesigns.

The Dutch roll mode (lateral, order 2) is basically a combined yaw androll oscillation. The period length and damping varies depending on theaircraft design. In civil aviation a weakly damped Dutch roll mode mustbe corrected by a yaw-damper.

The remaining 4 states can be considered to be pure integrators, as theaircraft dynamics are independent of its position on the earth and thedirection it is flying to. The altitude of flight, i. e., −z_(I), hasgot an influence on the air density and therefore on all aerodynamicforces. However, the dynamics are significantly slower than thementioned aircraft dynamics. Thus, it can be considered as a slowlyvarying parameter.

The following reduction of the system with 12 states to a kinetic modelwith 6 states is based on the following assumptions of a reducedtranslational kinetic model. Neglecting the dynamics of how the forcesthrust T, drag D, lift L and the roll angle φ originate. This must bejustified by using inner loop controllers which are considered to besufficiently faster than the outer loop. The thrust T acts perfectlyopposite to the drag D, i. e., the effective force in flight directionis T−D. In reality, depending on the angle of attack and the mounting ofthe propulsion unit, the effective direction of thrust force may not bealigned exactly opposite to the drag. The flight state is coordinated,i. e., the side force δF is zero. Hence, to obtain lateral forces in theinertial frame a rotation by φ is required which results in a lateralcomponent of lift L. Only steady wind is considered, i. e., for the windvelocity holds v_(W)≈const. Highly dynamic changes of wind in turbulentair are considered to act as external disturbance forces and torques.The inertial velocity v_(I) can be expressed as the sum ofin-air-velocity v_(A) and wind velocity v_(I)=v_(A)+v_(W) and

$\frac{{dv}_{I}}{dt} = {\frac{{dv}_{A}}{dt}.}$Henceforth, the index A indicates a representation with reference to theuniformly moving air mass. As the wind is considered to be constant, thetransformation from v_(I) to v_(A) is Galilean.

The velocity vector with reference to the uniformly moving air massv_(A) is described in the inertial frame using spherical coordinates

$\begin{matrix}{{v_{A} = \begin{bmatrix}{V_{A}\mspace{14mu}{\cos\left( \gamma_{A} \right)}\;{\cos\left( \psi_{A} \right)}} \\{V_{A}\mspace{14mu}{\cos\left( \gamma_{A} \right)}{\sin\left( \psi_{A} \right)}} \\{{- V_{A}}\mspace{14mu}{\sin\left( \gamma_{A} \right)}}\end{bmatrix}},} & (7)\end{matrix}$where V_(A) denotes the Euclidean norm ∥v_(A)∥₂, γ_(A) the climb angleand ψ_(A) the in-air flight direction of the aircraft. With thegravitational acceleration g, the derivatives can be found as

$\begin{matrix}{\frac{{dV}_{A}}{dt} = {\frac{T - D}{m} - {g\;{\sin\left( \gamma_{A} \right)}}}} & (8) \\{\frac{d\;\gamma_{A}}{dt} = \frac{{L\;{\cos(\varphi)}} - {m\; g\;{\cos\left( \gamma_{A} \right)}}}{m\; V_{A}}} & (9) \\{\frac{d\;\psi_{A}}{dt} = \frac{L\;{\sin(\varphi)}}{m\; V_{A}{\cos\left( \gamma_{A} \right)}}} & (10)\end{matrix}$

To obtain the derivative of the inertial position r_(I), i. e., thevelocity with reference to the inertial frame v_(I), the steady windvelocity v_(W)≈const. is taken into account, resulting in

$\begin{matrix}{\frac{{dr}_{1}}{dt} = {v_{I} = {v_{A} + v_{W}}}} & (11) \\{\frac{d^{2}r_{I}}{{dt}^{2}} = {\frac{{dV}_{A}}{dt} = {{\frac{\partial v_{A}}{\partial V_{A}}\frac{{dV}_{A}}{dt}} + {\frac{\partial v_{A}}{\partial\gamma_{A}}\frac{d\;\gamma_{A}}{dt}} + {\frac{\partial v_{A}}{\partial\psi_{A}}{\frac{d\;\psi_{A}}{dt}.}}}}} & (12)\end{matrix}$

Inserting (8, 9, and 10) into (12) one obtains the followingrepresentation:

$\begin{matrix}{\frac{d^{2}r_{I}}{{dt}^{2}} = {\begin{bmatrix}0 \\0 \\g\end{bmatrix} + {{D\left( {\gamma_{A},\psi_{A}} \right)}{\frac{1}{m}\begin{bmatrix}{T - D} \\{L\;{\sin(\varphi)}} \\{L\;{\cos(\varphi)}}\end{bmatrix}}}}} & (13) \\{{D\left( {\gamma_{A},\psi_{A}} \right)} = {\begin{bmatrix}{{\cos\left( \gamma_{A} \right)}{\cos\left( \psi_{A} \right)}} & {- {\sin\left( \psi_{A} \right)}} & {{- {\sin\left( \gamma_{A} \right)}}{\cos\left( \psi_{A} \right)}} \\{{\cos\left( \gamma_{A} \right)}{\cos\left( \psi_{A} \right)}} & {\cos\left( \psi_{A} \right)} & {{\sin\left( \gamma_{A} \right)}{\sin\left( \psi_{A} \right)}} \\{- {\sin\left( \gamma_{A} \right)}} & 0 & {- {\cos\left( \gamma_{A} \right)}}\end{bmatrix}.}} & (14)\end{matrix}$

The calculated state space model consists of 6 translational states,r_(I) and v_(I). The dynamics depend on the force inputs and theirorientation due to φ, γ_(A), and ψ_(A). γ_(A), and ψ_(A) describe theorientation of v_(A)=v_(I)−v_(W), they implicitly depend on v_(I) andv_(W). D being a rotation matrix of ψ_(A) and θ_(A), D is an orthogonalmatrix, implying det(D)=1. Therefore, D is regular for arbitrary ψ_(A)and γ_(A). As a first result, the aspired reduced translational kineticmodel is found as (13 and 14).

Based on the assumption that due to the inner loop controllers T, L, andφ can be manipulated independently, and m, g, D, γ_(A), and ψ_(A) areknown, the following transformation is made

$\begin{matrix}{\begin{bmatrix}\overset{\sim}{T} \\{\overset{\sim}{L}}_{S} \\{\overset{\sim}{L}}_{C}\end{bmatrix} = {\begin{bmatrix}{T - D} \\{L\;{\sin(\varphi)}} \\{L\;{\cos(\varphi)}}\end{bmatrix} = {m\;{D\left( {\gamma_{A},\psi_{A}} \right)}^{- 1}\left( {\begin{bmatrix}0 \\0 \\{- g}\end{bmatrix} + u} \right)}}} & (15)\end{matrix}$

This results in an exact linear behavior of the plant and reduces (13and 14) to

$\begin{matrix}{\frac{d^{2}r_{I}}{{dt}^{2}} = {u.}} & (16)\end{matrix}$

The actual desired values T, L, and φ to be realized by the innercontrol loops can be calculated as

$\begin{matrix}{T = {\overset{\sim}{T} + D}} & (17) \\{L = \sqrt{{\overset{\sim}{L}}_{C}^{2} + {\overset{\sim}{L}}_{S}^{2}}} & (18) \\{\varphi = {\arctan\left( \frac{{\overset{\sim}{L}}_{S}}{{\overset{\sim}{L}}_{C}} \right)}} & (19)\end{matrix}$

Hence, for transforming desired Cartesian accelerations u into [T Lφ]^(T), it is important to have accurate estimations of D, γ_(A), ψ_(A),and m.

The assumption of zero side force, corresponding to a suitably specifiedyaw angle, and freely definable roll angle implicates neglecting theDutch roll mode, the roll mode, and the spiral mode. The assumption offreely definable L, corresponding to a suitably specified pitch angle,implicates neglecting the short-period mode. Hence, in total 4 lateralstates and 2 longitudinal states are reduced to obtain a puretranslational kinetic model. Thus, the reduction in states can beinterpreted as neglecting the angular orientation of the aircraft, whichdetermines the aerodynamic forces, and only consider the states oftranslation.

Being an interchange of potential energy, correlated to the state z_(I),and kinetic energy, correlated to the velocity v_(I), the states of thephugoid mode are still part of the reduced kinetic model (13 and 14).This can be seen by specifying the lift L in (8, 9, 10) to beL=c_(L0)V_(A) ², c_(LO)=const., and thrust T=D, i. e., T compensates fordissipative drag D, resulting in

$\begin{matrix}{\frac{{dV}_{A}}{dt} = {{- g}\;{\sin\left( \gamma_{A} \right)}}} & (20) \\{\frac{d\;\gamma_{A}}{dt} = \frac{{c_{L\; 0}v_{A}^{2}\;{\cos(\varphi)}} - {m\; g\;{\cos\left( \gamma_{A} \right)}}}{m\; V_{A}}} & (21) \\{\frac{d\;\psi_{A}}{dt} = \frac{c_{L\; 0}v_{A}^{2}\;{\sin(\varphi)}}{m\; V_{A}{\cos\left( \gamma_{A} \right)}}} & (22)\end{matrix}$

For initial straight and level flight, i. e., φ=0 and the initial valuesγ_(A0)=0 and arbitrary ψ_(A0), the equilibrium airspeed results asV_(A0) according to L=mg=c_(L0)V_(A0) ². With ΔV_(A)=V_(A)−V_(A0),Δγ_(A)=γ_(A)−γ_(A0), and Δψ_(A)=ψ_(A)−ψ_(A0), linearizing the dynamics(8, 9, 10) results in

$\begin{matrix}{\frac{d\;\Delta\; V_{A}}{dt} = {{- g}\;{\Delta\gamma}_{A}}} & (23) \\{\frac{d\;{\Delta\gamma}_{A}}{dt} = {\frac{2_{g}}{V_{A\; 0}^{2}}\Delta\; V_{A}}} & (24) \\{\frac{d\;{\Delta\psi}_{A}}{dt} = 0.} & (25)\end{matrix}$

This is the characteristic form of an oscillator. For having compensatedthe only dissipative force D, the oscillation is undamped. With g=9.81m/s², the oscillation period, is

$\begin{matrix}{{T_{p\; h}\left( V_{A\; 0} \right)} = {{\frac{\sqrt{2}\pi}{g}V_{A\; 0}} = {0.453\frac{s^{2}}{m}{V_{A\; 0}.}}}} & (26)\end{matrix}$

The presented reduced kinetic model (13 and 14) can be of interest forairliners, where angle of attack α and side slip angle β are measuredand therefore γ_(A) and ψ_(A) can be calculated for known orientation ofthe aircraft. Furthermore, in many cases thrust T of the propulsionsystem can be determined and accurate models for estimating lift L anddrag D exist. Without the knowledge of these quantities, (15) and (17,18, 19) can't be calculated.

For not having available an accurate angle of attack and side slip anglemeasurement in low-cost UAVs, like the fixed wing aircraft used in thiswork for the validation of the PFC-concept, and most general aviationaircraft, the concept can't be used in this manner. Thus, the kineticmodel (13 and 14) is adapted hereinafter, pursuing the idea ofneglecting the dynamics of orientation to obtain a purely translationalkinematic model.

In the following a Reduced Translational Kinematic Model is introduced.If angle of attack and side slip angle measurements, necessary for thecalculation of γ_(A) and Vψ_(A), and knowledge of T, D, and L are notavailable, a further simplification can be done by considering only theresulting accelerations with reference to the body frame a_(B)=[a_(xB)a_(yB) a_(zB)]^(T) to be freely definable. Thus, the main difference tothe above is firstly the use of desired accelerations instead of desiredforces as inputs to the system. Secondly, the accelerations are orientedalong the body axes instead of the air flow axes. For example, the liftL is defined to be orthogonal to the undisturbed air flow pointing up,while a_(zB) is defined orthogonal to the longitudinal body axispointing down. In general, the two directions differ, e. g., due to theangle of attack.

The derivation is based on the following assumptions. Neglecting thedynamics of how the accelerations a_(xB), a_(zB) and the roll angle φoriginate. This must be justified by using inner loop controllers, thatare considered to be sufficiently faster than the outer loop. The flightstate is coordinated, i. e., a_(yB)=0. Hence, to obtain lateral forcesin the inertial frame a rotation by φ is required which results in alateral component of a_(zB). Only steady wind v_(w)≈const. isconsidered. a_(B) complies with accelerations measured by the utilizedaccelerometers which are not capable of measuring gravity, as discussed.Thus, the total acceleration acting on the aircraft is a_(I)=g+R_(I)^(B)a_(B), with g=[0 0 g]^(T).

By splitting R_(I) ^(B) into D(θ,ψ) and a canonical rotation about thebody x-axis by φ, the kinematic model results as

$\begin{matrix}{\frac{d^{2}r_{I}}{{dt}^{2}} = {\begin{bmatrix}0 \\0 \\g\end{bmatrix} + {{D\left( {\theta,\psi} \right)}\begin{bmatrix}a_{xB} \\{a_{zB}{\sin(\varphi)}} \\{a_{zB}{\cos(\varphi)}}\end{bmatrix}}}} & (27) \\{{D\left( {\theta,\psi} \right)} = {\begin{bmatrix}{{\cos(\theta)}{\cos(\psi)}} & {\sin(\psi)} & {{\sin(\theta)}{\cos(\psi)}} \\{{\cos(\theta)}{\sin(\psi)}} & {- {\cos(\psi)}} & {{\sin(\theta)}{\sin(\psi)}} \\{- {\sin(\theta)}} & 0 & {\cos(\theta)}\end{bmatrix}.}} & (28)\end{matrix}$

The hereby found kinematic model will also be the chosen representationfor the design of the path following control and the flight tests withthe fixed wing UAV. For D being a rotation matrix of yaw rotation ψ andthereafter pitch rotation θ, D is an orthogonal matrix, implyingdet(D)=1. Therefore, D is regular for arbitrary yaw and pitch angles.

Based on the assumption that due to the inner loop controllers a_(xB),a_(zB), and φ can be manipulated independently, the followingtransformation is made

$\begin{matrix}{\begin{bmatrix}a_{xB} \\a_{{zB},S} \\a_{{zB},C}\end{bmatrix} = {\begin{bmatrix}a_{xB} \\{a_{zB}{\sin(\varphi)}} \\{a_{zB}{\cos(\varphi)}}\end{bmatrix} = {{D\left( {\theta,\psi} \right)}^{- 1}{\left( {\begin{bmatrix}0 \\0 \\{- g}\end{bmatrix} + u} \right).}}}} & (29)\end{matrix}$

This results in an exact linear behavior of the plant and reduces (27and 28) to

$\begin{matrix}{\frac{d^{2}r_{I}}{{dt}^{2}} = u} & (30)\end{matrix}$

The actual desired values a_(zB) and φ to be realized by the innercontrol loops can be calculated as

$\begin{matrix}{{- a_{zB}} = \sqrt{a_{{zB},C}^{2} + a_{{zB},S}^{2}}} & (31) \\{\varphi = {\arctan\left( \frac{a_{{zB},S}}{a_{{zB},C}} \right)}} & (32)\end{matrix}$

To maintain the analogy to the lift L, mostly, −a_(zB) will be usedinstead of a_(zB). −a_(zB) is claimed to be positive, which means, thatdesired accelerations below weightlessness aren't considered in thiswork. φ is claimed to fulfill

${{\varphi } < \frac{\pi}{2}},$which means that inverted flight isn't considered either. Therefore, thesolution of (31 and 32) is unambiguous.

Finally, with (27 and 28) and the transformations (29), (31) and (32), arepresentation of a fixed wing aircraft is found, which perfectly suitsfor a path following control on the output r_(I) by means of inner loopcontrollers for a_(xB), a_(zB), φ, and a_(yB), which shall be zero for acoordinated flight.

The following aims to design inner loop controllers to reduce thedynamics of a fixed wing aircraft to the fundamental kinematic ones. Thebasic idea is to use directly measurable quantities like angular ratesand accelerations to obtain high bandwidth reference tracking anddisturbance rejection. This approach proves to be effective whendisturbances are hard to model or are unknown as it is the case, e. g.,for sticking friction. In this context, by measuring acceleration anddesigning a feedback control on the actuation force, the mentionednonlinear dynamics can be drastically simplified. The desired valuesφ^(des), θ^(des), a_(zB) ^(des), a_(yB) ^(des), V_(A) ^(des) from theouter loop are transformed to the control inputs δ_(A), δ_(E),δ_(F),δ_(R), δ_(T). The measured quantities V_(A), ω_(xB), ω_(yB), a_(zB),a_(yB), a_(xB) are another input for the corresponding controller, butare not indicated explicitly. The roll angle φ and pitch angle θ are notmeasured directly but estimated by the extended Kalman filter.

1. a) Identify the transfer function from δ_(A) to ω_(xB),

${i.e.},{P_{\omega_{xB}} = \frac{{\hat{w}}_{xB}}{{\hat{\delta}}_{A}}}$

-   -   b) Design a controller C_(ω) _(xB) to stabilize C_(ω) _(xB)        P_(ω) _(xB) by closing the feedback loop with the measurement        ω_(xB).    -   c) Design a controller C_(φ), to stabilize C_(φ)P_(φ), with

$P_{\varphi} = \frac{\hat{\varphi}}{{\hat{\omega}}_{xB}^{des}.}$2. a) Identify the transfer function from δ_(E) to ω_(yB),

${i.e.},{P_{\omega_{yB}} = \frac{{\hat{w}}_{yB}}{{\hat{\delta}}_{E}}}$

-   -   b) Design a controller C_(ω) _(yB) to stabilize C_(ω) _(yB)        P_(ω) _(yB) by closing the feedback loop with the measurement        ω_(yB).    -   c) Design a controller C_(θ) to stabilize C_(θ)P_(θ), with

$P_{\theta} = {\frac{\hat{\theta}}{{\hat{\omega}}_{yB}^{des}}.}$3. a) Identify the transfer function from δ_(F) to a_(zB),

${i.e.},{P_{a_{zB}} = {\frac{{\hat{a}}_{zB}}{{\hat{\delta}}_{F}}.}}$

-   -   b) Design a controller C_(a) _(zB) to stabilize C_(a) _(zB)        P_(a) _(zB) by closing the feedback loop with the measurement        a_(zB).        4. a) Identify the transfer function from δ_(R) to a_(yB),

${i.e.},{P_{a_{yB}} = {\frac{{\hat{a}}_{yB}}{{\hat{\delta}}_{R}}.}}$

-   -   b) Design a controller C_(a) _(yB) to stabilize C_(a) _(yB)        P_(a) _(yB) by closing the feedback loop with the measurement        a_(yB).        5. a) Identify the transfer function from δ_(T) to a_(xB),

${i.e.},{P_{a_{xB}} = {\frac{{\hat{a}}_{xB}}{{\hat{\delta}}_{T}}.}}$

-   -   b) Design a controller C_(a) _(xB) to stabilize C_(a) _(xB)        P_(a) _(xB) by closing the feedback loop with the measurement        a_(xB).    -   c) Design a controller C_(V) _(A) , to stabilize C_(V) _(A)        P_(V) _(A) , with

$P_{V_{A}} = \frac{{\hat{V}}_{A}}{{\hat{a}}_{xB}^{des}}$

In the above, the Laplace variable s is omitted for clarity, e.g., P_(ω)_(xB) represents P_(ω) _(xB) (s). The Laplace transform of quantities,which are also discussed in the time domain, is indicated by a hatsymbol above the variable, e.g., {circumflex over (ω)}_(xB). Thevariable and additional index z, e.g., P_(z,ω) _(xB) , indicates thetime discrete domain. Analogously, the variable and the additional indexq, e.g., P_(q,ω) _(xB) , indicates the Tustin domain.

The system, which is a solid body in the 3D-space with 6 degrees offreedom with 5 independent control inputs, is almost fully actuated.Especially by also using the flaps as dynamic input rather than only fortake-off and landing, the longitudinal subsystem is fully actuated,i.e., the 3 longitudinal degrees of freedom x_(I), z_(I), θ have thesame number of independent inputs δ_(E), δ_(F), δ_(T). Consequently, thelateral subsystem with the degrees of freedom y_(I), φ, ψ is missing oneinput to be fully actuated. δ_(A) primarily generates a roll moment tocontrol ω_(xB) and in further consequence φ. δ_(R) can be either used tocontrol ω_(zB) and in further consequence the yaw angle ψ, or to controla_(yB) by influencing the side slip angle β. During cruise flight,coordinated flight is achieved by a_(yB)=0. On the contrary, for landingthe alignment along the runway center line with orientation ψ_(RW) isimportant, i. e., ψ=ψ_(RW). As landing is not considered in this work,δ_(R) will be used to maintain coordinated flight, thus, will be used tocontrol a_(yB). Therefore, ω_(zB) is specified by the desired values ofa_(xB), a_(yB), a_(zB), ω_(xB), ω_(yB), i. e., cannot be controlledindependently, as the system is underactuated.

To obtain a fully actuated system, another input which directlygenerates side force would be necessary, e.g., kind of a wing invertical direction equipped with flaps. Consequently, a similar controlconcept as it will be demonstrated for the longitudinal system couldalso be implemented for the lateral system. A direct side force controland direct lift control can be used to control forces independently fromorientation and to improve the dynamic response of the forces. Fordirect lift control by means of flap actuation, the actuation time isthe limiting factor how fast lift can be changed. For state of the artairliners the flaps usually are only used for long-term modification ofthe lift characteristics, thus, extend slowly. The model aircraft whichwill be used for validation of the concept is equipped with the sameservos for flap actuation as for the other control surfaces. Therefore,highly dynamic flap extension is possible.

An adapted structure is investigated which allows to almost double theperformance of the a_(zB) control loop. Instead of the completelydecoupled SISO cascades, a dynamic filter output of δ_(F) is added toδ_(E). Furthermore, δ_(E) isn't used to realize ω_(yB) ^(des) anymorebut to return δ_(F) to its neutral position, e.g., δ_(F) ^(des)=0.Therefore, the inner control loop, where a_(xB), a_(yB), a_(zB), ω_(zB),and ω_(yB), are controlled, changes. ω_(yB) isn't controllableindependently anymore but specified by δ_(F) ^(des)=0. Thus, in thefinal concept, ω_(yB) and ω_(zB) aren't controlled but specified bya_(xB) ^(des), a_(yB) ^(des)=0, a_(zB) ^(des), ω_(xB) ^(des), and δ_(F)^(des)=0.

The identification task is accomplished by the application of a chirpsignal on the respective input. By means of the identification of thetransfer function P_(ω) _(yB) from the input δ_(E) to the output ω_(yB),the necessary calculations are explained. As it is rarely the case thatthe air is completely smooth, multiple identification patterns are flownand averaged to reduce the effect of disturbances. Still, a reasonablysmooth air is desirable to obtain satisfying results. Another importantconsideration is the frequency band to be identified which shouldcontain the dynamics of interest. For practical reasons, the time of onepattern is limited to approximately 20 s to keep the aircraft closeenough to the remote pilot to maintain visual contact. Therefore, thefrequency band shouldn't be chosen too wide, as this would result inweak excitation of the corresponding frequencies. On the other hand,choosing the frequency band too narrow results in only littleinformation about the frequency dependency of the dynamics. It turnedout to be a good trade-off to choose one decade in frequency around thearea of interest. Another important consideration is the amplitude ofthe signal. If the amplitude is too large, the state of equilibrium onwhich the linear identification is based may be left excessively. On theother hand, a very low amplitude doesn't excite the system enough toachieve differentiation from disturbance and noise.

Another important precondition for the identification process is thechosen equilibrium state on which the chirp signal is superposed. Thisinformation is specified for each identification procedure. In case ofthe identification of P_(ω) _(yB) the point of equilibrium and the taskof the other controllers during the identification is given in FIG. 2.As ω_(yB) represents a longitudinal state, lateral dynamics aresuppressed by maintaining coordinated, straight flight, i. e., a_(yB)=0m/s² and wings leveled, i. e., φ=0 rad. As the control action of theV_(A)-control could have parasitic effects on the identification, duringthe identification pattern throttle is kept constant, i. e.,δ_(T)=δ_(T,0)=const. The value δ_(T,0) is determined by an onlinelow-pass filter of first order of δ_(T), which determines the mean valuebefore the start. For ω_(yB)-control, the flaps aren't used and kept at0. The identification input δ_(E,ID) consists of the chirp signalitself, a trim value δ_(E,0), determined analogously to δ_(T,0), and thepossibility of additional manual input. The controllers for φ, θ,a_(yB), and V_(A) can be the result of a preceding controller design ormanually tuned PI-controllers.

The response is analyzed in the discrete time domain by assuming atransfer function of third order. This order may be different for thetransfer functions of the other identification procedures. In the caseof P_(ω) _(yB) it is expected to be identified in the region of theshort-period mode, an oscillation of second order. Additionally, thebandwidth of reference tracking of the servos is limited to similarfrequencies, which is taken into account by another state. Thus, thestructure of the time-discrete plant can be written as

$\begin{matrix}{{P_{z,\omega_{yB}}(z)} = \frac{{b_{1}z^{- 1}} + {b_{2}z^{- 2}} + {b_{3}z^{- 3}}}{1 + {a_{1}z^{- 1}} + {a_{2}z^{- 2}} + {a_{3}z^{- 3}}}} & (33)\end{matrix}$

With T_(s)=0.02 s, y_(k)=ω_(yB)(kT_(s)), and u_(k)=δ_(E)(kT_(s)), in thediscrete time domain, P_(z,ω) _(yB) corresponds toy _(k) =b ₁ u _(k−1) +b ₂ u _(k−2) +b ₃ u _(k−3) −a ₁ y _(k−1) −a ₂ y_(k−2) −a ₃ y _(k−3)  (34)

Based on measured values y₁, y₂, . . . , y_(N−1) and u₁, u₂, . . . ,u_(N−1), summarized in the data matrix

$\begin{matrix}{s = \begin{bmatrix}{- y_{3}} & {- y_{2}} & {- y_{1}} & u_{3} & u_{2} & u_{1} \\{- y_{4}} & {- y_{3}} & {- y_{2}} & u_{4} & u_{3} & u_{2} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{- y_{N - 1}} & {- y_{N - 2}} & {- y_{N - 3}} & u_{N - 1} & u_{N - 2} & u_{N - 3}\end{bmatrix}} & (35)\end{matrix}$and with y=[y₄ y₅ . . . y_(N)]^(T), an optimal solution for theparameter vector p=[a₁ a₂ a₃ b₁ b₂ b₃]^(T), according to the leastsquares problem

$\begin{matrix}{{\min\limits_{p}{\left( {y - {y_{est}(p)}} \right)^{T}\left( {y - {y_{est}(p)}} \right)}},{y_{est} = {Sp}}} & (36)\end{matrix}$is found to bep*=(S ^(T) S)⁻¹ S ^(T) y.  (37)

Finally, with p*=[0.41-0.90 1.01-1.66 1.12-0.33]^(T) the estimated plantP_(z,ω) _(yB) can be written as

$\begin{matrix}{P_{z,\omega_{yB}} = {\frac{{{- 1.66}z^{- 1}} + {1.12z^{- 2}} - {0.33z^{- 3}}}{1 + {0.41z^{- 1}} - {0.90z^{- 2}} + {1.01z^{- 3}}}.}} & (38)\end{matrix}$

In the following the procedure how the linear controllers are designedis shown. As an example, the cascaded pitch-controller to track areference input θ^(des) for the pitch angle is considered. The innerloop controller C_(ω) _(yB) generates the desired value of elevatordeflection δ_(E), which is realized by a servomotor, to track thedesired rate ω_(yB) ^(des) of the higher-level controller C_(θ).

The design of C_(ω) _(yB) is based on the beforehand identified plantP_(z,ω) _(yB) . By applying the bilinear transformation

$\begin{matrix}{{P_{q,\omega_{yB}}(q)} = {\frac{{{- 0.5571}q^{3}} + {84.66q^{2}} - {4150q} - {1.255 \cdot 10^{5}}}{q^{3} + {59.29q^{2}} + {2774q} + {4.567 \cdot 10^{4}}}.}} & (39)\end{matrix}$to (38), the representation of the system in the Tustin domain isdetermined as

$z = \frac{1 + {{qT}_{s}/2}}{1 - {{qT}_{s}/2}}$

For this plant a PI-controller is designed by using two tuningparameters ω_(0dB) and ω_(P) in the form

$\begin{matrix}{{C_{q,\omega_{yB}}(q)} = {{V\left( \omega_{0d\; B} \right)}{\frac{1 + \frac{q}{\omega\; p}}{q}.}}} & (40)\end{matrix}$

The parameter ω_(0dB) determines the frequency where the open looptransfer function L_(q,ω) _(yB) =C_(q,ω) _(yB) P_(q,ω) _(yB) is supposedto cross the 0 dB level. The gain factor V(ω_(0dB)) is adjusted toachieve this property. ω_(0dB) roughly characterizes the bandwidth.ω_(p) specifies the root of the nominator of C_(q,ω) _(yB) (q) andtherefore the frequency, where integral-dominant action changes toproportional-dominant action, i.e., where a phase-shift of −90° changesto 0°.

A good trade-off between bandwidth and robustness is found for ω_(0dB)=9rad/s and ω_(p)=3ω_(0dB)=27 rad/s. Thus, at ω_(0dB) the integral actionof the PI-controller dominates, which is necessary to obtain a slope of−20 dB per decade due to the almost constant magnitude of the plantP_(z,ω) _(yB) itself. The −20 dB slope is a design criterion which isrequired to avoid creeping settling of the controlled quantity to thedesired reference value. The gain factor V(ω_(0dB)) follows as

$\begin{matrix}{{{L_{q,{\omega_{y}B}}\left( {j\;\omega_{0d\; B}} \right)}} = {{{V\left( \omega_{0d\; B} \right)}{{\frac{1 + \frac{j\;\omega_{0d\; B}}{\omega_{p}}}{j\;\omega_{0d\; B}}{P_{q,{\omega_{y}B}}\left( {j\;\omega_{0d\; B}} \right)}}}}\overset{!}{=}1}} & (41) \\{{{V\left( \omega_{0d\; B} \right)}{{\frac{1 + \frac{j\;\omega_{0d\; B}}{\omega_{p}}}{j\;\omega_{0d\; B}}{P_{q,{\omega_{y}B}}\left( {j\;\omega_{0d\; B}} \right)}}}^{- 1}} = {3.265.}} & (42)\end{matrix}$

Thus, the controller in the Tustin domain results as

$\begin{matrix}{{C_{q,{\omega_{y}B}}(q)} = {3.265{\frac{1 + {q/\left( {27\frac{rad}{s}} \right)}}{q}.}}} & (43)\end{matrix}$

At ω_(0dB) the magnitude is confirmed to be 1 while the phase plot showsa phase reserve ϕ_(0dB)=180°−119.530=60.470. By applying the inversebilinear transformation

${q = {\frac{2}{T_{s}}\frac{z - 1}{z + 1}}},$the corresponding representation in the discrete time domain is found tobe

$\begin{matrix}{{C_{z,\omega_{yB}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.1536},{k_{I} = 3.2648}} & (44)\end{matrix}$

For a simulated closed loop step response for ω_(yB) ^(des)=1 rad/s, thecontroller effort of elevator deflection δ_(E) shows a steady statevalue of 0.364 and a maximum value of 0.4. Thus, assuming linearbehavior, the restrictions for the elevator |δ_(E)|≤1 are fulfilled for

${{{\omega_{yB}^{des}} \leq \frac{1}{0.4}} = 2.5},$which will be chosen to be the output restriction for the higher-levelcontroller. Furthermore, the gain factor of C_(ω) _(yB) is corrected bythe factor

$\frac{V_{ref}}{V_{A}}$according to V_(ref)=12 m/s. The transfer function from {circumflex over(δ)}_(E) to {circumflex over (ω)}_(yB) shows an amplification which isdirectly proportional to the airspeed V_(A). This known relation istherefore compensated by the division with V_(A). The same correction isalso done for the roll rate ω_(xB) for the roll-controller cascade.

For the design of the higher-level controller C_(θ), the plant

$P_{\theta} = \frac{\hat{\theta}}{{\hat{\omega}}_{yB}^{des}}$must be determined. For this purpose, one option is to carry outidentification flights with the loop of C_(ω) _(yB) closed, ω_(yB)^(des) as input, and θ as response for the identification.

Assuming that the performance of the controller C_(ω) _(yB) correspondsto the designed one, another possibility exists. For φ=0 rad, ω_(yB) isthe derivative of θ. For higher bank angles φ≠0 rad, this relationdoesn't hold any longer and would have to be corrected. However, thepitch controller is only used for identification maneuvers in straightflight conditions, thus the assumption

$\omega_{yB} = \frac{d\;\theta}{dt}$is justified. Due to

${\omega_{yB} = \frac{d\;\theta}{dt}},{i.e.},{\hat{\theta} = \frac{{\hat{\omega}}_{yB}}{s}}$and the assumption that

${T_{\omega_{yB}} = {\frac{{\hat{\omega}}_{yB}}{{\hat{\omega}}_{yB}^{des}} = \frac{c_{\omega_{yB}}P_{\omega_{yB}}}{1 + {c_{\omega_{yB}}P_{\omega_{yB}}}}}},$as designed, the transfer function from {circumflex over (ω)}_(yB)^(des) to {circumflex over (θ)} can be calculated to be

${P_{\theta} = {\frac{\hat{\theta}}{{\hat{\omega}}_{yB}^{des}} = {{\frac{\hat{\theta}}{{\hat{\omega}}_{yB}}\frac{{\hat{\omega}}_{yB}}{{\hat{\omega}}_{yB}^{des}}} = {\frac{1}{s}T_{\omega_{yB}}}}}},{i.e.},$it can directly be determined by the complementary sensitivity functionof the control design of C_(ω) _(yB) . Analogous considerations can bemade for the cascaded controllers C_(φ), C_(ω) _(xB) , and C_(V) _(A) ,C_(a) _(xB) .

Analogously to the PI-controller design of C_(ω) _(yB) , with ω_(0dB)=3rad/s and ω_(p)=0.2ω_(0dB)=0.6 rad/s, the controller C_(θ) results as

$\begin{matrix}{{C_{z,\theta} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 2.9570},{k_{I} = {1.7636.}}} & (45)\end{matrix}$

In analogous manner to the identification procedure and design for thepitch-controller cascade, the roll-controller cascade consisting ofC_(ω) _(xB) and C_(φ) is designed. First of all, the plant P_(ω) _(xB)is identified. The chosen equilibrium state and controller conditionsbefore and during the identification pattern can be found in FIG. 3.δ_(A) is driven according to the identification pattern and theadditional possibility to make manual corrections to keep level flight.δ_(E) is controlled to keep a constant pitch angle θ^(des)=0 rad tomaintain the equilibrium state. Otherwise, especially for increasingbank angles, the aircraft's nose would drop, which leads to airspeedincrease and a spiral turn. δ_(R), δ_(T) are kept constant during theidentification pattern to avoid parasitic effects.

The identified transfer function P_(z,ω) _(xB) is

$\begin{matrix}{P_{z,\omega_{xB}} = {\frac{{{- 0.24}z^{- 1}} + {0.75z^{- 2}} - {0.17z^{- 3}}}{1 - {1.85z^{- 1}} + {1.25z^{- 2}} - {0.30z^{- 3}}}.}} & (46)\end{matrix}$

Based on this transfer function with ω_(0dB)=10 rad/s andω_(P)=3ω_(0dB)=30 rad/s, a PI-controller C_(z,ω) _(yB) is designed inthe form

$\begin{matrix}{{C_{z,\omega_{xB}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.1278},{k_{I} = {2.9485.}}} & (47)\end{matrix}$

For a simulated closed loop step response for ω_(xB) ^(des)=1 rad/s, thecontroller effort of aileron deflection δ_(A) shows a steady state valueof 0.31 and a maximum value of 0.36. Thus, assuming linear behavior, therestrictions for the aileron |δ_(A)|≤1 are fulfilled for |ω_(xB)^(des)|≤1/0.36=2.78, which will be chosen to be the output restrictionfor the higher-level controller.

The plant P

$P_{\varphi} = \frac{\hat{\varphi}}{{\hat{\omega}}_{xB}^{des}}$is found in an analogous manner as P_(θ). Assuming that the performanceof the controller C_(ω) _(xB) is according to the designed one, and dueto

${\omega_{xB} = \frac{d\;\varphi}{dt}},{i.e.},{\hat{\varphi} = \frac{{\hat{\omega}}_{xB}}{s}},$the transfer function from ω_(xB) ^(des) to φ can be calculated as

$P_{\varphi} = {\frac{\hat{\varphi}}{{\hat{\omega}}_{xB}^{des}} = {{\frac{\hat{\varphi}}{{\hat{\omega}}_{xB}^{des}}\frac{{\hat{\omega}}_{xB}}{{\hat{\omega}}_{xB}^{des}}} = {\frac{1}{s}{T_{\omega_{xB}}.}}}}$

Based on this transfer function with ω_(0dB)=3 rad/s andω_(p)=0.2ω_(0dB)=0.6 rad/s, a PI-controller C_(φ) is designed in theform

$\begin{matrix}{{C_{z,\varphi} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 2.9442},{k_{I} = {1.7560.}}} & (48)\end{matrix}$

The transfer function from {circumflex over (δ)}_(A) to {circumflex over(ω)}_(xB) shows an amplification which is directly proportional to theairspeed V_(A). This known relation is therefore compensated by thedivision with V_(A) in the form

$\frac{V_{ref}}{V_{A}}.$

An a_(yB)-controller is designed for the plant P_(a) _(yB) with therudder δ_(R) as input. FIG. 4 shows the conditions for theidentification. The identified transfer function P_(z,a) _(yB) is

$\begin{matrix}{P_{z,a_{yB}} = \frac{{{- 0.43}z^{- 1}} - {1.12z^{- 2}} + {1.83z^{- 3}}}{1 - {1.21z^{- 1}} - {0.035z^{- 2}} + {0.28z^{- 3}}}} & (49)\end{matrix}$

Based on this transfer function with ω_(0dB)=2 rad/s andω_(p)=3ω_(0dB)=6 rad/s, a PI-controller C_(a) _(yB) is designed

$\begin{matrix}{{C_{z,a_{yB}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.0415},{k_{I} = {0.2347.}}} & (50)\end{matrix}$

The transfer function from {circumflex over (δ)}_(R) to â_(yB) shows aquadratic dependency on the airspeed V_(A). This known relation istherefore compensated by the division with V_(A) ² in the form

$\frac{V_{ref}^{2}}{V_{A}^{2}}.$

In analogous manner to the identification procedure and design for thepitch-controller cascade, the airspeed-controller cascade consisting ofC_(a) _(xB) and C_(V) _(A) is designed. The chosen equilibrium state andcontroller conditions before and during the identification pattern canbe found in FIG. 16. The plant P_(a) _(xB) is identified as describedabove.

The identified transfer function P_(z,a) _(xB) is

$\begin{matrix}{P_{z,a_{xB}} = {\frac{{0.054z^{- 1}} + {0.44z^{- 2}} - {0.85z^{- 3}} + {0.76z^{- 4}}}{1 - {1.23z^{- 1}} + {0.41z^{- 2}} - {0.16z^{- 3}} + {0.044z^{- 4}}}.}} & (51)\end{matrix}$

Based on this transfer function with ω_(0dB)=10 rad/s andω_(P)=3ω_(0dB)=30 rad/s, a PI-controller C_(a) _(xB) is designed

$\begin{matrix}{{C_{z,a_{xB}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.1817},{k_{I} = {0.8651.}}} & (52)\end{matrix}$

The plant

$C_{V_{A}} = \frac{{\hat{V}}_{A}}{{\hat{a}}_{xB}^{des}}$is found in an analogous manner as P_(θ). Assuming that the performanceof the controller C_(a) _(xB) is according to the designed one and dueto

${a_{xB} = \frac{{dV}_{A}}{dt}},{i.e.},{{\hat{V}}_{A} = \frac{{\hat{a}}_{xB}}{s}},$the transfer function from a_(xB) ^(des) to V_(A) can be calculated as

$P_{V_{A}} = {\frac{{\hat{V}}_{A}}{{\hat{a}}_{xB}^{des}} = {{\frac{{\hat{V}}_{A}}{{\hat{a}}_{xB}}\frac{{\hat{a}}_{xB}}{{\hat{a}}_{xB}^{des}}} = {\frac{1}{s}{T_{a_{xB}}.}}}}$

Based on this transfer function with ω_(0dB)=1 rad/s andω_(P)=0.2ω_(0dB)=0.2 rad/s, a PI-controller C_(V) _(A) is designed

$\begin{matrix}{{C_{z,V_{A}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.9860},{k_{I} = {0.1968.}}} & (53)\end{matrix}$

In the following, a method is developed to control a_(zB) by means ofactuating the wing, i. e., δ_(F), not only actuates the flaps but alsolowers the ailerons. Thus, the curvature of the whole span is changed toinfluence the aerodynamic accelerations generated by the wing. Furtherdevelopments could improve this technique by using adaptive materials tochange the wing shape more precisely than by means of 3 discrete controlsurfaces.

Generally, the control design could be done in the same manner asdemonstrated above. However, a pitch oscillation of the short-periodmode results in an antiresonance of vertical acceleration −a_(zB) whenexcited by δ_(F). If the ω_(yB) control was fast enough, it wouldsuppress this oscillation. Still, as the performance of the systemlimits the achievable bandwidth of the ω_(yB) control, a differentmethod is developed to overcome the antiresonance and therefore improvethe performance.

To obtain constant magnification in the control loop, instead ofcontrolling a_(zB) itself the coefficient

$\begin{matrix}{c_{zB} = {\frac{a_{zB}}{V_{A}^{2}}\frac{V_{ref}^{2}}{a_{ref}}}} & (54)\end{matrix}$is calculated, i. e., the quadratic dependency on the airspeed V_(A) iscompensated. Due to scaling with a_(ref)=g=9.81 m/s² and V_(ref)=12 m/s,at level flight conditions with V_(A)=V_(ref) follows −c_(zB)=1. As forthis control concept the indicated airspeed V_(A) is expected to bealways higher than a minimum airspeed V_(A)>V_(A,min)>0 m/s, thedivision is well defined. The desired independence of c_(zB) on V_(A) isconfirmed by identification procedures at different airspeeds. To thisend, a series of test flights is executed where the response of c_(zB)to the chirp of δ_(F) indeed doesn't alter significantly for differentairspeeds.

To understand the difference in lift generation due to δ_(E) and δ_(F)deflection, the model of lift force (5 and 6) is repeated

$\begin{matrix}{L = {c_{L}\frac{\rho}{2}V_{A}^{2}S}} & (55) \\{{c_{L \approx}c_{L\; 0}} + {c_{L,\alpha}\alpha} + {c_{L,q}q} + {c_{L,\delta_{E}}\delta_{E}} + {c_{L,\delta_{F}}{\delta_{F}.}}} & (56)\end{matrix}$

The differences between c_(L) and the introduced c_(zB) are theorientation and a constant scaling factor. Lift L is defined to beorthogonal to the undisturbed airflow and positive in the upwardsdirection. The body acceleration a_(zB) is orthogonal to thelongitudinal body axis and positive to the downwards direction. Thedifferent scaling is influenced by the mass m, for L being a force anda_(zB) an acceleration, other factors like wing area S, and theintroduced normalization by a_(ref) and V_(ref). For small angles ofattack α, the fundamental characteristic for lift generation isconsidered to be similar for c_(L) and c_(zB), thus, also c_(zB) can bewritten asc _(zB) ≈c _(zB,0) +c _(zB,α) α+c _(zB,q) q+c _(zB,δ) _(E) δ_(E) +c_(zB,δ) _(F) δ_(F).  (57)

Reference is made to FIG. 5. 202 shows |P_(q,δ) _(E) |, plot 204 shows|P_(q,δ) _(F) |, plot 206 shows arg(P_(q,δ) _(F) ) and plot 208 showsarg(P_(q,δ) _(E) ). By identifying

$P_{\delta_{E}} = {- \frac{{\hat{c}}_{zB}}{{\hat{\delta}}_{E}}}$with δ_(F)=0, and

$P_{\delta_{F}} = {- \frac{{\hat{c}}_{zB}}{{\hat{\delta}}_{F}}}$with δ_(E)=δ_(E,0)=const., the main difference of the two methods togenerate lift gets obvious.

Generation of lift by the elevator δ_(E) is mainly based on changing theangle of attack α, i. e., the effect of c_(zB,α)α dominates. Thecomponent c_(zB,δ) _(E) δ_(E) is a parasitic effect of the elevator dueto the force at the horizontal tail, which is necessary to generate thepitching moment to change the angle of attack. In the case of atailplane, this parasitic force results to be opposite to the resultinglift of the wing. As the dynamics of the force at the tailplane onlydepend on the actuator of the elevator, while the angle of attack mainlychanges according to the short-period mode, a nonminimum phasecharacteristic can be observed. Therefore, generation of lift by theelevator inevitably entails a lag, which depends on the short-periodmode and in case of a tailplane the polarity of the force at first evenresults to be opposite to the desired one. The simulated step responsefor the identified transfer function of −c_(zB) excited by δ_(E)illustrates this characteristic in FIG. 6, plot 210 shows the stepresponse of −c_(zB) due to step input δ_(E) and plot 212 shows the stepresponse of −c_(zB) due to step input of δ_(F). FIG. 5 emphasizes thelow-pass characteristic of the short-period mode, with a characteristicfrequency of ω_(SP)=13.15 rad/s.

Generation of lift by flap deflection δ_(F) directly influences thecomponent c_(zB,δ) _(F) δ_(F) of (57). If no parasitic pitching momentis caused by the wing deformation, the stationary value of the angle ofattack α remains the same, i. e., c_(zB,α)α doesn't change for lowexcitation frequencies. Still a dynamic effect is observed, which causesan a oscillation. This oscillation causes an antiresonance approximatelyin the region of ω_(SP) where the short-period mode has got a phaseshift of 90°. As illustrated by FIG. 7 this phase shifted a oscillationresults in an effect c_(zB,α)α which cancels out part of the desiredcomponent c_(zB,δ) _(F) δ_(F). Namely, if δ_(F) by means of c_(zB,δ)_(F) δ_(F) excites a sinusoidal acceleration a_(zB,δ) _(F) , asinusoidal trajectory 214 will set in, where maximum upwardsacceleration is found at minimum height at point B in FIG. 7. Thedirection parallel to the trajectory is indicated by e_(∥). Consideringthe orientation of the aircraft, for low frequencies the pitch anglefollows the change of orientation of the trajectory. This is due to thecharacteristic of the aircraft to orient into the airflow, i. e., tokeep alpha constant. For higher frequencies the aircraft orientationdoesn't follow the trajectory direction ideally anymore. According tothe decoupling of the short-period mode, the longitudinal axis which isindicated by e_(xB) oscillates with a phase shift, i. e., theorientation of e_(xB) lags behind the orientation of e_(∥). FIG. 7illustrates this lag of the aircraft's longitudinal direction e_(xB)behind the trajectory direction e_(∥) for a 90° phase shift, whichoccurs at the resonance frequency of the short-period mode. Therefore,e_(xB) reaches the ideal orientation into the airflow of point A, i. e.,the orientation of e_(∥,A), as late as at point B. This causes anegative angle of attack at point B and therefore a downwardsacceleration a_(zB,α) which cancels out part of the flap effect a_(zB,δ)_(F) . The vertical movement in FIG. 7 is overstated to make the effectclear. Also, the rotation of a_(zB) due to the rotated body x-axise_(xB) isn't illustrated, without affecting the underlying principle.This effect is also confirmed by simulation of a reduced model with thestates θ, ω_(yB), V_(A), γ_(A). Thereby the longitudinal dynamics, whichresult in the short-period mode and the phugoid mode can beinvestigated. x_(I) and z_(I) have no influence on the dynamics and canbe neglected like the lateral states. As expected, the suppressioneffect of the antiresonance depends on the damping of the short-periodmode which determines the resulting magnitude of the α oscillation inthe region of the resonance frequency. The identified plant

$- \frac{{\hat{c}}_{zB}}{{\hat{\delta}}_{F}}$shows the described antiresonance in FIG. 5.

FIG. 6 shows the step response of −c_(zB) to a step in δ_(F) of theidentified plant P_(δ) _(F) , where higher frequency dynamics have beenremoved. Therefore, low pass behavior of the sensors and actuators isneglected. This is done to emphasize the potential of immediate responseof vertical acceleration by flap deflection. In contrast, the responseof −c_(zB) to a step in δ_(E) at first shows nonminimal characteristicand due to the short-period mode needs about 0.4 s to rise. This delayin acceleration generation by δ_(E) can't be improved by fasteractuators, as the aircraft in any case needs to rotate to alter theangle of attack. However, by improving the actuator dynamics of δ_(F) aswell as the system set-up concerning time delay, a very fast response ofvertical acceleration to δ_(F) can be achieved. For the givenexperimental set-up the advantage is limited by the overall time delayof the system. Still, more than double the bandwidth is possiblecompared to c_(zB) control by δ_(E), as long as the antiresonance iscompensated.

To achieve this, a dynamic filter is designed to compensate theoscillation by use of δ_(E). To this end, the input δ_(F) is filteredand the output of the dynamic filter is added to δ_(E). Thus, theintention is to reduce the α oscillation by means of anticipating theeffect of δ_(F) and compensate it using δ_(E). Hence, a dynamic filterF_(δ) _(E) _(,δ) _(F) is designed to compensate the describedantiresonance by means of adding a feedforward elevator deflection{circumflex over (δ)}_(E,FF)=F_(δ) _(E) _(,δ) _(F) {circumflex over(δ)}_(F) to the control input {circumflex over (δ)}_(E). Thus, assuminglinear superposition, the new resulting transfer function {tilde over(P)}_(δ) _(F) from {circumflex over (δ)}_(F) to −ĉ_(zB) becomes

$\begin{matrix}{{- {\hat{c}}_{zB}} = {{{P_{\delta_{F}}{\hat{\delta}}_{F}} + {P_{\delta_{E}}{\hat{\delta}}_{E,{FF}}}} = {\underset{\underset{\mspace{121mu}{\overset{\sim}{P}}_{\delta_{F}}\mspace{101mu}}{︸}}{\left( {P_{\delta_{F}} + {P_{\delta_{E}}F_{\delta_{E},\delta_{F}}}} \right)}\mspace{11mu}{\hat{\delta}}_{F}}}} & (58)\end{matrix}$

As ansatz for the filter F_(δ) _(E) _(,δ) _(F) , a linear filter ofsecond order is chosen. As the filter shall only operate at higherfrequencies, where the antiresonance appears, the steady state responseshall be 0. Thus, the ansatz in the Tustin-domain is

$\begin{matrix}{{F_{q,{\delta\; E},{\delta\; F}}(q)} = {\frac{{b_{1}q} + {b_{2}q^{2}}}{q^{2} + {2\;\xi_{F}\omega_{F}q} + \omega_{F}^{2}}.}} & (59)\end{matrix}$

Using the Matlab command fminunc( ), an optimization problem is solvedto find an optimal fit of {tilde over (P)}_(δ) _(F) to a desired plant{tilde over (P)}_(δ) _(F) ^(des) by tuning the parameter vector [b₁b₂ξ_(F) ω_(F)]. The cost function to be minimized is chosen to be

$\begin{matrix}{J = {\sum\limits_{\omega_{k}{\epsilon\Omega}}{{{{\overset{\sim}{P}}_{q,\delta_{F}}\left( {j\;\omega_{K}} \right)} - {{\overset{\sim}{P}}_{q,\delta_{F}}^{des}\left( {j\;\omega_{k}} \right)}}}}} & (60)\end{matrix}$

With Ω=2π{0.3; 0.4; . . . ; 3.9; 4}. Thus, the magnitude response isoptimized at 38 equally spaced frequencies in the considered frequencyband. The phase response isn't included in the cost function, as thetime delay of the system wasn't determined explicitly and can't becompensated by a causal filter. For a desired plant {tilde over (P)}_(δ)_(F) ^(des)=1.45, the optimized filter parameters result as b₁=18.15,b₂=0.968, ξ_(F)=0.991, ω_(F)=20.04 rad/s. FIG. 8 depicts the Bode plotsof the original plant |P_(q,δ) _(F) | 222 and arg(P_(q,δ) _(F) ) 230,the desired plant |{tilde over (P)}_(q,δ) _(F) ^(des)| 216 andarg({tilde over (P)}_(q,δ) _(F) ^(des)) 228, the filter |P_(q,δ) _(E)_(,δ) _(F) | 224 and arg(F_(q,δ) _(E) _(,δ) _(F) ) 226, and theresulting plant |P_(q,δ) _(F) | 220 and arg(P_(q,δ) _(F) ) 232.

For implementation, the filter designed in the Tustin-domain istransformed to the discrete domain resulting in

$\begin{matrix}{F_{z,\delta_{E},\delta_{F}} = {\frac{0.7997 - {1.347\; z^{- 1}} + {0.5472\; z^{- 2}}}{1 - {1.336z^{- 1}} + {0.4474z^{- 2}}}.}} & (61)\end{matrix}$

FIG. 10 illustrates the effect of the filter by comparing the responsesof −c_(zB) 236, 240 to the identification pattern for δ_(F) 234, 238without and with the filter F_(z,δ) _(E) _(,δ) _(F) , which weremeasured during test flights. In the time interval from 2 s until 12 sfrequencies in the region of 2 Hz are excited. This corresponds to thefrequency band where the antiresonance occurs, cf. FIG. 8. Theantiresonance is clearly evident in the response without F_(δ) _(E)_(,δ) _(F) , which is illustrated in the upper plot of FIG. 10. Thelower plot of FIG. 10 shows the response of −c_(zB) including the filterF_(δ) _(E) _(,δ) _(F) . The antiresonance is confirmed to becompensated, i. e., almost constant magnification can be achieved.

The plant {tilde over (P)}_(δ) _(F) is identified, with identificationconditions of the table of FIG. 9, as

$\begin{matrix}{{\overset{\sim}{P}}_{z,\delta_{F}} = {\frac{0.137 + {0.135z^{- 1}}}{1 - {1.25z^{- 1}} + {0.44z^{- 2}}}.}} & (62)\end{matrix}$

For this plant the c_(zB) controller is designed, as demonstrated forthe controller C_(ω) _(yB) . With ω_(0dB)=15 rad/s and ω_(P)=1.8ω_(0dB)=27 rad/s, the controller C_(z,c) _(zB) results as

$\begin{matrix}{{C_{z,c_{zB}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = 0.4773},{k_{I} = {10.1479.}}} & (63)\end{matrix}$

The effective range of accelerations that can be generated by δ_(F) islimited. FIG. 11 shows the closed loop step response of −c_(zB) 242, andδ_(F) 244, for −c_(zB) ^(des)=1, which corresponds to−a_(zB)=a_(ref)=9.81 m/s² at V_(A)=V_(ref)=12 m/s. The flaps 244 show asteady state deflection of δ_(F)=0.7. Furthermore, flaps parasiticallygenerate drag for high deflections. Thus, a method is implemented toreturn δ_(F) to δ_(F) ^(des)=0. FIG. 12 shows the structure 300 of theimplemented control to track a desired c_(zB) ^(des). δ_(E) consists ofthe sum created by an adder 304 of the filter output of F_(δ) _(E) _(,δ)_(F) of a feed forward filter 302, which is designed to compensate theantiresonance, and the output (candidate elevator angle δ_(E) ^(C)) of afirst controller 306, C_(δ) _(F) , which will be designed to countersteady state flap deflection. C_(δ) _(F) will be designed ofconsiderably lower bandwidth than a second controller 308, C_(c) _(zB) ,to prevent interactions of the two controllers 306, 308. For slow flightδ_(F) ^(des) can also be specified to δ_(F) ^(des)>0, however, the inputlimits must be considered to maintain enough control input range ofδ_(F) for C_(c) _(zB) .

Thus, the response P_(δ) _(F) _(,δ) _(E) of δ_(E) to δ_(F) isidentified, while C_(c) _(zB) is tracking c_(zB) ^(des)=c_(zB,0), asindicated in table of FIG. 13. c_(zB,0) is calculated by a low-passfilter analogously to δ_(T,0). The transfer function is identified as

$\begin{matrix}{P_{z,{\delta\; F},{\delta\; E}} = \frac{0.124 - {0.197\; z^{- 1}}}{1 - {1.61z^{- 1}} + {0.634z^{- 2}}}} & (64)\end{matrix}$

For ω_(0dB)=3 rad/s, a fifth of the bandwidth of C_(c) _(zB) , andω_(P)=2ω_(0dB)=6 rad/s, the second controller 308 C_(δ) _(F) results as

$\begin{matrix}{{C_{z,\delta_{F}} = {k_{P} + {k_{I}\frac{T_{s}}{z - 1}}}},{k_{P} = {- 0.1948}},{k_{I} = {- 1.1025}}} & (65)\end{matrix}$

As C_(c) _(zB) produces an output disturbance that shall be returned tothe desired value δ_(F) ^(des)=0, the response to a disturbance ofδ_(F)=1 is presented in FIG. 14 instead of reference tracking, whereinFIG. 14 shows a response of δ_(E) 246 and δ_(F) 248 to an outputdisturbance δ_(F)=1. An overshoot of about 6% and a rise time of about0.7 s can be observed.

Finally, the functionality of the control concept is illustrated in FIG.15. The diagram shows a typical lift curve 252, which mainly generates−c_(zB), over the angle of attack α. The effect of the flaps δ_(F) 250,254 is to shift the curve −c_(zB)(α), upwards for positive δ_(F). Theeffect of the elevator δ_(E) is to manipulate α to reach differentpoints of −c_(zB)(α), which constitutes the conventional method forvarying the lift. The introduced control concept includes the flaps forhighly dynamic c_(zB) variations. FIG. 15 outlines the setpoint changefrom −c_(zB,0) 256 to −c_(zB,1) 258. First, C_(c) _(zB) via δ_(F) with adesigned bandwidth of about 15 rad/s increases −c_(zB) in a highlydynamic manner. Then C_(δ) _(F) via δ_(E) returns δ_(F) to the initialvalue, in this case δ_(F) ^(des)=0. The different dynamics areemphasized by the small arrows, which shall correspond to constant timesteps. Thus, the new value −c_(zB,1) is reached with the high dynamicsof C_(c) _(zB) , whereas the horizontal movement approximately on thelevel of −c_(zB,1) is based on a fifth of the bandwidth, thus, 3 rad/s.The transition from δ_(F) to δ_(E) is similar to FIG. 14.

Not including δ_(F), thus, controlling c_(zB) directly via δ_(E), cf.P_(q,δ) _(E) in FIG. 5, would limit the bandwidth to approximately 7rad/s. Due to the low-pass characteristic of the short-period mode aboveapproximately 12 rad/s even faster actuators and smaller time delayscould not improve a direct control concept with δ_(E). This is due tothe fact that the whole aircraft must rotate to alter the lift by meansof changing α.

In contrast, the implementation of the introduced two-stage controlconcept proves to not be limited by the short-period mode. Therefore, afaster system could further rise the achievable bandwidth.

The curve −c_(zB)(α) is pretty linear until the region around α_(max)where the wings begin to stall and lift decreases again. By setting alimitation to a value below −c_(zB,max), e. g., −c_(zB,1), a simple, buteffective stall prevention can be implemented.

As a summary, various linear controllers are designed to fulfill twobasic tasks. Firstly, the pitch θ, roll φ, airspeed V_(A), and sideacceleration a_(yB) controllers are used to maintain coordinated,straight, and level flight at a desired airspeed. This is used formaintaining an equilibrium state during identification maneuvers andfurthermore, achieving a well defined initial state at the beginning oftest flights for validation of the PFC concept. Secondly, when the pathfollowing control is active, the inner loop controllers realize thedesired input u of the reduced plant

$\frac{d^{2}r_{I}}{{dt}^{2}} = {u.}$The results of (29), (31), (32), and (54) are combined to obtain thenecessary transformations to calculate a_(xB) ^(des), c_(zB) ^(des) andφ^(des), for given θ, ψ, g, and the desired input u as

$\begin{matrix}{\begin{bmatrix}a_{xB}^{des} \\a_{{zB},S}^{des} \\a_{{zB},C}^{des}\end{bmatrix} = {{D\left( {\theta,\psi} \right)}^{- 1}\left( {\begin{bmatrix}0 \\0 \\{- g}\end{bmatrix} + u} \right)}} & (66) \\{{- a_{zB}^{des}} = \sqrt{\left( a_{{zB},C}^{des} \right)^{2} + \left( a_{{zB},S}^{des} \right)^{2}}} & (67) \\{c_{zB}^{des} = {\frac{a_{zB}^{des}}{V_{A}^{2}}\frac{V_{ref}^{2}}{a_{ref}}}} & (68) \\{\varphi^{des} = {{atan}\left( \frac{a_{{zB},S}^{des}}{a_{{zB},C}^{des}} \right)}} & (69)\end{matrix}$with D according to (28). Thus, the tasks of C₇, C_(c) _(zB) , C_(a)_(xB) are to track the according desired value, which result in theinertial accelerations

$\frac{d^{2}r_{I}}{{dt}^{2}} = {u.}$The desired reference input for C_(δ) _(F) is chosen as δ_(F) ^(des)=0,which may be adapted during take-off and landing. The desired referenceinput for C_(a) _(yB) is a_(yB) ^(des)=0, according to coordinatedflight. By considering the influence of airspeed variation, the designedcontrollers may be used for the entire speed range. A common systemvariation in aviation is varying mass due to changing payload. Thenature of the acceleration controllers for a_(xB) and c_(zB), whichcorresponds to a_(zB), is the generation of forces. According to theconservation of linear momentum for constant mass, the forces areproportional to the accelerations by multiplication of mass. Thus, massvariations can be corrected by adapting the open loop gain of thecorresponding controller by the factor

$\frac{m}{m_{ref}},$where m_(ref) is the mass for which the controller has been designed.The actual mass m can be measured or approximately calculated during thestandard mass and balance calculations during flight preparation. As aconclusion, the inner loop controllers are designed to obtain a system

$\frac{d^{2}r_{I}}{{dt}^{2}} = u$by means of the nonlinear transformation (66, 67, 68, 69). Thenonlinearities of the aircraft are primarily based on cosine and sinefunctions of orientation and the dependencies of forces and torques onairspeed. These relations are known and can be compensated to obtain anexactly linear plant by means of inner loop control.

Based on the system

${\frac{d^{2}r_{I}}{{dt}^{2}} = u},$cf. (16) and (30), a path following control for a fixed wing aircraft isdesigned. Thus, by the assumed inner loop controllers, which realize thedesired accelerations a_(I) ^(des)≈u, this system appears in Brunovskyform with vector relative degree {2, 2, 2} for the output r_(I)=[x_(I)y_(I) z_(I)]^(T). The so-called path parameter ζ defines the desiredpath σ(ζ) in parameterized form. σ(ζ) is required to be twofolddifferentiable, and, if the resulting accelerations shall be continuous,twofold continuously differentiable. Henceforth, the dependency of timer_(I)(t), ζ(t), etc. is omitted for clarity, except for specialemphasis. The derivatives are abbreviated as

$\sigma^{\prime} = \frac{\partial\sigma}{\partial\zeta}$and {dot over (ζ)}=dζ/dt, etc.

The control objectives of PFC are defined as follows.

-   -   The output r_(I) converges asymptotically to the path σ(·),

${i.e.},\left. {\inf\limits_{\zeta}{{r_{I} - {\sigma(\zeta)}}}}\rightarrow 0 \right.$for t→∞. (asymptotic convergence)

-   -   If at t₀ the aircraft position r_(I) is on the path and the        velocity v_(I) is parallel to the path, i. e., ∃ζ₀:        r_(I)=σ(ζ₀)∧v_(I)=κσ′(ζ₀),κ∈        , then

${\inf\limits_{\zeta}{{{r_{I}(t)} - {\sigma(\zeta)}}}} = 0$for all t≥t₀. (invariance property)

-   -   The requirements for ζ(t) depend on the operation mode.        (tangential motion)

In the following, the operation mode for realizing a path followingcontroller with a desired path speed is deduced. Considering the patherror e_(p)=r_(I)−σ(ζ), the error dynamics are

$\begin{matrix}{e_{P} = {r_{I} - \underset{\underset{r_{P}}{︸}}{\sigma(\zeta)}}} & (70) \\{{\frac{d}{dt}e_{P}} = {{{\frac{d}{dt}r_{I}} - {\frac{d}{dt}\left( {\sigma(\zeta)} \right)}} = {v_{I} - {\underset{\underset{v_{P}}{︸}}{\sigma^{\prime}(\zeta)}\overset{.}{\zeta}}}}} & (71) \\{{\frac{d^{2}}{{dt}^{2}}e_{P}} = {{{\frac{d^{2}}{{dt}^{2}}r_{I}} - {\frac{d^{2}}{{dt}^{2}}\left( {\sigma(\zeta)} \right)}} = {u - {\underset{a_{P}}{\underset{︸}{\left( {{{\sigma^{''}(\zeta)}{\overset{.}{\zeta}}^{2}} + {{\sigma^{\prime}(\zeta)}\overset{¨}{\zeta}}} \right)}}.}}}} & (72)\end{matrix}$

For the special case on hand, also analogies to a PID-Control can bemade. Therefore, the coefficients of the feedback terms are named k_(P),k_(I), and k_(D). For

$\begin{matrix}{{u = {{{\sigma^{''}(\zeta)}{\overset{.}{\zeta}}^{2}} + {{\sigma^{\prime}(\zeta)}\overset{¨}{\zeta}} - {k_{P}e_{P}} - {k_{D}\frac{d}{dt}e_{P}}}},} & (73)\end{matrix}$the desired error dynamics

$\begin{matrix}{{{\frac{d^{2}}{{dt}^{2}}e_{P}} + {k_{D}\frac{d}{dt}e_{P}} + {k_{P}e_{P}}} = 0} & (74)\end{matrix}$

can be specified by tuning the parameters k_(P)>0 and k_(D)>0 accordingto the characteristic polynomial p(s)=s²+k_(D)s+k_(P). It is alsopossible to obtain individual dynamics for the different Cartesian axesby using the matrix K_(p)=diag([k_(p1) k_(p2) k_(p3)]) instead of thescalar k_(P), and K_(D) in an analogous manner. Furthermore, individualdynamics for other than the Cartesian axes of the inertial frame couldbe defined by use of nondiagonal K_(p) and K_(D). In this disclosure,just the case of scalar parameters is considered.

It seems reasonable to compare the designed error dynamics (74) to athree dimensional spring-damper-system, where the equilibrium pointcontinuously changes according to the path parameter ζ(t) and thecorresponding path position σ(ζ(t)). Thus, the terms in (73) can beinterpreted as spring action −k_(P)e_(p), damper action

${k_{D}\frac{d}{dt}e_{P}},$acceleration feedforward σ″(ζ){dot over (ζ)}² due to path parameterspeed {dot over (ζ)}, and acceleration feedforward σ′(ζ){umlaut over(ζ)} due to path parameter acceleration {umlaut over (ζ)}.

By introducing an integral error

${{\frac{d}{dt}e_{I}} = e_{P}},$the final control law for the input results as

$\begin{matrix}{u = {{{\sigma^{''}(\zeta)}{\overset{.}{\zeta}}^{2}} + {{\sigma^{\prime}(\zeta)}\overset{¨}{\zeta}} - {k_{P}e_{P}} - {k_{D}\frac{d}{dt}e_{P}} - {k_{I}e_{I}}}} & (75)\end{matrix}$

Therefore, the error dynamics extend to

$\begin{matrix}{{{{\frac{d^{3}}{{dt}^{3}}e_{P}} + {k_{D}\frac{d^{2}}{{dt}^{2}}e_{P}} + {k_{P}\frac{d}{dt}e_{P}} + {k_{I}e_{P}}} = 0},} & (76)\end{matrix}$with the additional tuning parameter k_(I)>0 yielding the correspondingcharacteristic polynomial p(s)=s³+k_(D)s²+k_(P)s+k_(I).

For the calculation of the control law (75), the course of the pathparameter ζ(t) must be specified, which depends on the operation mode.The operation mode is characterized by a freely definable path speedV_(P) ^(des)(t). In the general case of a path without naturalparameterization, the path parameter derivatives {dot over (ζ)} and{umlaut over (ζ)} have to be adapted continuously to achieve the desiredpath speed ∥V_(P)∥λ=V_(P) ^(des) (t), where λ ∈ {−1, 1} specifies theflying direction along the path, i. e., λ=sgn(V_(P) ^(des))=sgn({dotover (ζ)}). According to the desired path position r_(P)=σ(ζ), thedesired path speed follows as ∥v_(P)∥λ=∥σ′(ζ)∥{dot over (ζ)} and thus

$\begin{matrix}{\overset{.}{\zeta} = \frac{V_{P}^{des}(t)}{{\sigma^{\prime}(\zeta)}}} & (77) \\{v_{P} = {\frac{\sigma^{\prime}(\zeta)}{{\sigma^{\prime}(\zeta)}}{{V_{P}^{des}(t)}.}}} & (78)\end{matrix}$

For the thereby specified path parameter speed {dot over (ζ)} and withthe desired path acceleration A_(P) ^(des)(t)={dot over (V)}_(P)^(des)(t), the path parameter acceleration follows in the form

$\begin{matrix}{\overset{¨}{\zeta} = {\frac{A_{P}^{des}(t)}{{\sigma^{\prime}(\zeta)}} - {\frac{\left( {\sigma^{\prime}(\zeta)} \right)^{T}{\sigma^{''}(\zeta)}}{{{\sigma^{\prime}(\zeta)}}^{2}}{{\overset{.}{\zeta}}^{2}.}}}} & (79)\end{matrix}$

From these quantities the desired acceleration

$a_{P} = {\frac{d^{2}}{{dt}^{2}}r_{P}}$can be calculated asa _(P)=σ″(ζ){dot over (ζ)}²+σ′(ζ){umlaut over (ζ)}.  (80)If σ(ζ) is twofold continuously differentiable and V_(P) ^(des) (t) iscontinuously differentiable, i.e., σ″(ζ) and A_(P) ^(des)(t) arecontinuous, also a_(P) results to be continuous. For the case of anatural parameterization, the relations (σ′(ζ))^(T)σ′(ζ)=1 and therefore(σ′(ζ))^(T)σ″(ζ)=0 hold. Thus, {dot over (ζ)} and {umlaut over (ζ)} in(77) and (79) are reduced to{dot over (ζ)}=V _(P) ^(des)(t)  (81){umlaut over (ζ)}=A _(P) ^(des)(t).  (82)

Finally, by application of the path following control law (75), theinput for the inner controllers u, cf. (66, 67, 68, 69), is calculated.The operation mode to fly with a definable path speed V_(P) ^(des)(t) isachieved by choosing the path parameter ζ according to (77) and (79).

The invention claimed is:
 1. A method for controlling a pitch momentgenerator element of an aircraft, the method comprising: a step ofdetermining a pitch control input δ_(E) for the at least one pitchmoment generator element, based on a lift control input δ_(F) for atleast one lift generator element; wherein the step of determining thepitch control input δ_(E) based on the lift control input δ_(F) includesdetermining an output S_(E,FF)=F_(q,δE,δF)(q)δF of a feed forward filterF_(q,δE,δF)(q) for the lift control input δ_(F), and adding the outputδ_(E,FF) of the feed forward filter F_(q,δE,δF)(q) to a candidate pitchcontrol input δ_(E) ^(C), wherein the candidate pitch control inputδ_(E) ^(C) is controlled by one of a pilot, a flight controller and anautopilot; wherein determining the output δ_(E,FF) of the feed forwardfilter F_(q,δE,δF)(q) includes at least one of the following: high passfiltering the lift control input δ_(F) by a high pass; high passfiltering the lift control input δ_(F) by a lead-filter; high passfiltering the lift control input δ_(F) by a high pass of second order.2. The method according to claim 1, wherein the pitch control input forat least one pitch moment generator element is constituted by acommanded elevator angle δ_(E) to be determined based on a flap elementangle δ_(F) constituting a lift control input for at least one liftgenerator element.
 3. The method according to claim 1, wherein the liftcontrol input δ_(F) is commanded by one of a pilot, a flight controllerand an autopilot.
 4. The method according to claim 1, wherein the liftcontrol input δ_(F) for the at least one lift generator elementdetermines at least one of: a deflection of a flap; a deflection of aflaperon; a deflection of a spoiler; a deflection of a slat; an angle ofincidence of a wing; a lift generating change of a wing shape; a liftgenerating thrust of a micro-thruster; a lift generating thrust of apropulsion unit; and the pitch control input δ_(E) for the at least onepitch moment generator element determines at least one of: a deflectionof an elevator; an angle of incidence of a horizontal stabilizer; anangle of incidence of a canard foreplane; a pitch moment generatingchange of a wing shape; a pitch moment generating thrust of amicro-thruster; a pitch moment generating thrust of a propulsion unit.5. The method according to claim 1, wherein the feed forward filterF_(q,δE,δF)(q) is determined by the following formula:${F_{q,{\delta\; E},{\delta\; F}}(q)} = \frac{{b_{1}q} + {b_{2}q^{2}}}{q^{2} + {2\xi_{F}\omega_{F}q} + \omega_{F}^{2}}$wherein q is the Laplace variable; b₁ is a first optimization parameter;b₂ is a second optimization parameter; ξ_(F) is a third optimizationparameter; and ω_(F) is a fourth optimization parameter.
 6. The methodaccording to claim 1, further comprising the steps of: determining adesired lift control input δ_(F) ^(des); and determining the candidatepitch control input δ_(E) ^(C) based on the lift control input δ_(F) andthe desired lift control input δ_(F) ^(des) of the aircraft, wherein thestep of determining a candidate pitch control input δ_(E) ^(C) includesa PI control based on the lift control input δ_(F) and the desired liftcontrol input δ_(F) ^(des) of the aircraft.
 7. The method according toclaim 1, further comprising the steps of: determining an actual verticalacceleration a_(zB) of the aircraft; determining a desired verticalacceleration a_(zB) ^(des) of the aircraft; and controlling the liftcontrol input δ_(F) based on the actual vertical acceleration a_(zB) ofthe aircraft and the desired vertical acceleration a_(zB) ^(des) of theaircraft.
 8. The method according to claim 1, comprising pre-actuatingthe lift generator element thereby to subsequently allow it to bothdecrease lift and increase lift.
 9. A computer program product that whenloaded into a memory of a computer having a processor executes a methodfor controlling a pitch moment generator element of an aircraft, themethod comprising: a step of determining a pitch control input δ_(E) forthe at least one pitch moment generator element, based on a lift controlinput δ_(F) for at least one lift generator element; wherein the step ofdetermining the pitch control input δ_(E) based on the lift controlinput δ_(F) includes determining an output δ_(E,FF)=F_(q,δE,δF)(q)δF ofa feed forward filter F_(q,δE,δF)(q) for the lift control input δ_(F),and wherein determining the output δ_(E,FF) of the feed forward filterF_(q,δE,δF)(q) includes at least one of the following: high passfiltering the lift control input δ_(F) by a high pass; high passfiltering the lift control input δ_(F) by a lead-filter; and high passfiltering the lift control input δ_(F) by a high pass of second order.10. A flight controller adapted to control at least one pitch momentgenerator element of at least one of an aircraft and a component of aflight simulator, the flight controller being comprised in at least oneof an aircraft and a flight simulator, and comprising: a feed forwardfilter F_(q,δE,δF)(q) adapted to determine a pitch control input δ_(E)for the at least one pitch moment generator element, based on a liftcontrol input δ_(F) for the at least one lift generator element, bydetermining an output δ_(E,FF)=F_(δE,δF)(q)δF of the feed forward filterF_(q,δE,δF)(q) for the lift control input δ_(F), wherein the feedforward filter includes at least one of the following: a high passfilter for filtering the lift control input δ_(F); a lead filter forfiltering the lift control input δ_(F); a high pass of second order forfiltering the lift control input δ_(F).
 11. The flight controlleraccording to claim 10, wherein the pitch control input for at least onepitch moment generator element is constituted by a commanded elevatorangle δ_(E) to be determined based on a flap element angle δ_(F)constituting a lift control input for at least one lift generatorelement.
 12. The flight controller according to claim 10, wherein thelift control input is commanded by one of a pilot, a flight controllerand an autopilot.
 13. The flight controller according to claim 10,wherein the lift control input δ_(F) for the at least one lift generatorelement determines at least one of: a deflection of a flap; a deflectionof a flaperon; a deflection of a spoiler; a deflection of a slat; anangle of incidence of a wing; a lift generating change of a wing shape;a lift generating thrust of a micro-thruster; a lift generating thrustof a propulsion unit; and the pitch control input δ_(E) for the at leastone pitch moment generator element determines at least one of: adeflection of an elevator; an angle of incidence of a horizontalstabilizer; an angle of incidence of a canard foreplane; a pitch momentgenerating change of a wing shape; a pitch moment generating thrust of amicro-thruster; a pitch moment generating thrust of a propulsion unit.14. A flight controller according to claim 10, wherein the feed forwardfilter F_(q,δE,δF)(q) is determined by the following formula:${F_{q,{\delta\; E},{\delta\; F}}(q)} = \frac{{b_{1}q} + {b_{2}q^{2}}}{q^{2} + {2\xi_{F}\omega_{F}q} + \omega_{F}^{2}}$wherein q is the Laplace variable; b₁ is a first optimization parameter;b₂ is a second optimization parameter; ξ_(F) is a third optimizationparameter; and ω_(F) is a fourth optimization parameter.
 15. The flightcontroller according to claim 10, comprising an adder for adding anoutput δ_(E,FF) of the feed forward filter F_(q,δE,δF)(q) to a candidatepitch control input δ_(E) ^(C), wherein the candidate pitch controlinput δ_(E) ^(C) is controlled by one of a pilot, a flight controllerand an autopilot.
 16. The flight controller according to claim 15,further comprising a first controller adapted to determine a candidatepitch control input δ_(E) ^(C) based on the lift control input δ_(F) anda desired lift control input δ_(F) ^(des) of the aircraft, wherein thefirst controller is a PI controller.
 17. The flight controller accordingto claim 10, further comprising a second controller adapted to controlthe lift control input δ_(F) based on an actual vertical accelerationa_(zB) of the aircraft and a desired vertical acceleration a_(zB) ^(des)of the aircraft.
 18. The flight controller according to claim 10,comprising a configuration allowing it to pre-actuate the lift generatorelement thereby to subsequently allow it to both decrease lift andincrease lift.